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Engineering Traffic Networks for More Than One 
Busy Hour 


By M. EISENBERG 
(Manuscript received July 6, 1976) 


A procedure is described that is used to engineer traffic networks for 
more than one hour of point-to-point load data. The procedure differs 
significantly from existing methods, which are based upon the concept 
of “economic load on the last trunk” (ECCS). (When the peak-load 
hours on most routes coincide, however, the procedure reduces to the 
ECCS method.) This “multihour” procedure has been implemented 
in a computer program used in design studies of three end offices in the 
Los Angeles local network. For the cases examined, the multihour 
technique produced networks whose costs averaged approximately 7 
percent below those achieved with the presently used single-hour 
methods. Thus, the multihour technique appears to promise consid- 
erable cost benefits in future network designs. 


I. INTRODUCTION 


In this paper, we describe a procedure used to engineer networks for 
more than one hour of point-to-point traffic data. Specifically, for a given 
routing structure, set of switching and transmission costs, and point- 
to-point offered load between each pair of offices for each of several 
hours, this method produces a (nearly) least-cost network that satisfies 
the constraint that the blocking probabilities on all final groups be below 


1 


a predetermined value (the “grade of service”) for all hours. This mul- 
tihour procedure is a major departure from currently used single-hour 
methods based upon the concept of “economic load on the last trunk”! 
(ECCS). The new technique reduces to the ECCS method, however, 
when the peak-load hours on most routes coincide. After computer 
programs are written and operating practices are developed, this new 
procedure should become suitable for routine field use. 

The underlying theoretical basis for multihour engineering was de- 
veloped by Rapp.” Rather than attempting to construct an optimal so- 
lution, however, Rapp proposed an alternative approximate technique. 
Our aim is to get an exact solution. Although we do not fully achieve this 
aim, we obtain significant improvement in network performance relative 
to a single-hour approach. 

A computer program that implements the multihour procedure was 
used to study three end offices located in the Los Angeles area. For the 
cases examined, where a significant amount of noncoincidence of 
peak-load hours existed, the multihour method produced network cost 
savings averaging 7 percent over the single-hour methods currently 
employed. In addition, in each case, a very sizable reduction of tandem 
switching load was achieved. 


ll. SINGLE-BUSY-HOUR ENGINEERING 


Before discussing the multihour technique, let us first review the 
considerations involved in engineering for a single hour. Figure 1 depicts 
a single high-usage group, the direct route, overflowing to an alternate 
route. (For now, we make the simplifying assumption that the alternate 
route consists only of a single trunk group. We later consider more re- 
alistic alternate-route configurations.) The cost per trunk of the direct 
route is Cp, and the cost per trunk of the alternate route is C4 (which 
is assumed to include the cost of tandem switching). The offered load 
in the hour being considered is a. The problem is to determine the value 
of n, the number of trunks in the high-usage group, so that the total cost 
is minimized; however, the minimization of cost is subject to the con- 
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Fig. 1—High-usage trunk group overflowing to alternate route. 
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Fig. 2—Single-hour trunk-group sizing. 


straint that the blocking probability on the alternate route is below a 
predetermined value. 

The total cost is equal to the cost of trunks on the direct route plus the 
cost of trunks on the alternate route. To simplify matters, it can be as- 
sumed that the alternate-route cost is composed of a fixed component 
(the cost of carrying the “background” alternate-route load) plus a 
variable component (the cost of carrying the overflow) whose magnitude 
is proportional to the amount of overflow traffic.* Since the cost required 
to carry the background alternate-route load is independent of n, we may 
neglect this component and write the cost to be minimized as 


COST = CA Re) + Cpn. (1) 
Y 


Here, y is the marginal capacity of the alternate route, and B(n,a) is the 
Erlang-B blocking probability. The marginal capacity is the amount of 
additional traffic that it is assumed can be offered to the alternate route, 
at fixed blocking, for the addition of one trunk. Thus, aB(n,qa) is the load 
overflowing to the alternate route, 1/y aB(n,a) is the assumed number 
of additional alternate-route trunks required to carry this overflow, and 
C,4/y aB(n,a) is the cost of these additional trunks. Cpn is, of course, 
the cost of trunks on the direct route. These two components of cost and 
their sum are shown as a function of n in Fig. 2. The total cost is seen to 
be a U-shaped curve having a minimum at the point indicated. This point 
is determined by the condition that the rate of change of COST with 
respect to n be equal to zero: 


* This assumption is not quite true, particularly when the peakedness of the overflow 
traffic is taken into account. Nevertheless, in most cases of interest, the assumption yields 
a configuration whose cost differs negligibly from the optimal plan. 
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£ cost = 0. (2) 
dn 
From eq. (1), this implies 
d y 
== [aR =e 


where Cr = C4/Cp is the cost ratio, as discussed in Ref. 1. The quantity 
on the left-hand side of eq. (3), the rate of change of overflow with respect 
to the number of trunks on the high-usage group, is very nearly equal 
to the “load on the last trunk.” * The quantity on the right-hand side 
of eq. (3) is the “economic” load on the last trunk, or ECCS (CCS is 100 
call seconds per hour). Thus, the minimum cost is achieved by sizing the 
high-usage group such that its load on the last trunk is equal to its 
“economic” value, y/Cr. 

In this discussion, it is assumed that the network is designed to carry 
only a single hour’s load. In practice, of course, the load on the high-usage 
group, as well as the background load on the alternate route, varies from 
hour to hour. The question arises as to which of the hours of loads should 
be used to engineer the group. 

It is clear that it would be uneconomical to engineer a high-usage group 
for its individual group busy hour if this hour does not coincide with the 
busy hour of the alternate route; the alternate route has spare capacity 
in off-hours. A moment’s reflection reveals that the appropriate hour 
for which to size the group is the alternate-route’s busy hour; only in this 
hour does the cost of carrying the overflow traffic from the high-usage 
group need to be considered. 

This fact has long been recognized by traffic engineers.! The method 
of choosing the engineering hour which was adopted, consequently, in- 
volved the concept of the “cluster busy hour.” A “cluster” is defined as 
a set of high-usage trunk groups originating at the same office and 
overflowing to a common alternate-route leg, together with the alter- 
nate-route leg itself. The cluster busy hour is defined as that time-con- 
sistent hour for which the total load offered to the cluster (specifically, 
the sum of the carried loads on all high-usage groups in the cluster, plus 
the offered load on the alternate-route leg) is maximum. It was assumed 
that the alternate-route busy hour would be the same as the cluster busy 
hour, and thus the adopted engineering practice was to size every high- 
usage group for its cluster-busy-hour load. 

The difficulties that can arise with this method, however, are illus- 
trated by the example in Fig. 3. In the figure, we show a simple network 


* “Toad on the last trunk” is defined to be a[B(n-1,a)-B(n,a)| if the group has n 
trunks. 
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cluster consisting of two one-way outgoing high-usage groups, A and B, 
overflowing to a common alternate route, F. The cost ratio is assumed 
to be 2, and the marginal capacity of the alternate route is assumed to 
be 28 CCS. Since a total of 1100 CCS is offered to the cluster in Hour 1 
and 1000 CCS is offered in Hour 2, the cluster busy hour is Hour 1. For 
the Hour-1 loads the ECCS method yields a network consisting of 10 and 
26 trunks on high-usage groups A and B, respectively. It is seen, however, 
that this network has a high overflow from group A during Hour 2. This 
high overflow occurs because the group is engineered for only 300 CCS, 
while in Hour 2 the offered load is 600 CCS. As a consequence, the total 
load on the alternate route is greater during Hour 2, contradicting the 
original assumption that the alternate route was busier in Hour 1. To 
guarantee a given grade of service in both hours, it is necessary to add 
trunks to the alternate route for its load in Hour 2.* Under the as- 
sumption that the number of extra alternate-route trunks required to 
carry this load is nr = 276/y = 9.8, we find the total cost of the network 
to be $55.60.7 

Figure 3 also shows the network derived on the basis of the Hour-2 
loads. For this network, due to high overflow from group B during Hour 
1, the final-group busy hour is Hour 1, again contradicting the initial 
assumption. The total cost of this network is $59.40. 

The third network in Fig. 3 was derived using the multihour technique. 
As can be seen, this network nearly equalizes the load on the alternate 
route in the two hours. (In general, the multihour technique tends to 
equalize the hourly loads on the alternate route or routes.) The cost of 
this network is $46.80, substantially less than either single-hour net- 
work. 

This example illustrates some of the problems inherent in single-hour 
engineering methods and the potential improvement obtainable with 
the multihour technique. We describe this technique in detail in Section 
II. 


lil. MULTIHOUR ENGINEERING 


Figure 4 again shows a single high-usage group overflowing to an al- 
ternate route, where the trunk costs Cp and C4 are defined as before. 
The loads a, and a2 are offered to the high-usage group in Hours 1 and 
2, respectively. Also shown in the figure are the background loads in 


* In practice, the servicing up of the alternate route might take place after the engineered 
a was in operation, when the service degradation in the side hour was actually ob- 
served. 

t The use of an assumed marginal capacity for the alternate route, while reasonable for 
the purpose of sizing high-usage groups, is actually inappropriate for determining trunk 
requirements on the alternate route. Our aim here, however, is merely to obtain a rough 
indication of alternate-route cost for comparative purposes. 
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TOTAL COST = 10+ 26+2x9.8 
= 55.6 








TOTAL COST = 204+134+2x 13.2 
= 59.4 


=> np = 13.2 


TOTAL COST = 14+20+2x6.4 





a 165 => n= 6.4 


Fig. 3—Comparison of engineering methods in the presence of noncoincidence. 


Hours 1 and 2, A, and Ag, offered to the alternate route. The background 
loads are the total loads offered to the alternate route, not including the 


overflow from the high-usage group under consideration. 
In sizing the high-usage group, we attempt to minimize total cost. We 
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Fig. 4—High-usage group with offered load for two hours. 
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must recognize, however, that since the grade of service must be guar- 
anteed for both hours, the cost of the alternate route depends upon the 
greater of its total offered loads in Hours 1 and 2. The cost is thus given 
by the formula* 


Ay + a,B(n,a1) 


+ Cpn. 4 
Ao + aeB(n,ae) oe (4) 


cost = &4 Max | 
Y 


The total load on the alternate route during Hour 1 is equal to the 
background load A, plus the overflow from the high-usage group, 
a,B(n,a,). The total load during Hour 2, similarly, is Az + a2B(n,ay). 
The controlling load for the alternate route is the greater of these. The 
total cost equals the maximum alternate-route load times C4/y plus the 
cost of the high-usage group Cpn.* 

The term “multihour engineering” denotes the process of designing 
a network by searching along the cost curve of eq. (4) for each group—or 
actually, the more general cost curve of eq. (5) discussed below—to de- 
termine the minimum-cost point. The optimal number of high-usage 
trunks in each group determined by the use of this technique varies 
depending upon the loads and costs. Figure 5 shows two different cases 
that can arise. 

In Case I, the plot of cost curves shows that the Hour-1 load on the 
alternate route dominates the Hour-2 load for all n. In this case, there- 
fore, the maximization operator of eq. (4) always selects Hour 1, as 
suggested by the shading of this cost curve in the figure. Thus, in Case 
I, the multihour method reduces exactly to the single-hour method by 
using the Hour-1 load. This example illustrates the case discussed above, 
where the use of the cluster-busy-hour concept yields the correct solu- 
tion. Note that the correct solution is confirmed if the actual alternate- 
route busy hour, determined by examining the load offered after engi- 
neering, is the same as that originally assumed. 

Case II in Fig. 5 illustrates a different type of behavior. In this example, 
the background load on the alternate route is greater in Hour 2, and the 
offered load on the high-usage group is greater in Hour 1. Thus, for n 
small there is heavy overflow during Hour 1, and this causes the total 
alternate-route load to be greater in Hour 1. For n large, however, the 
overflow is small, so that the alternate-route load is greater in Hour 2 
due to the background component. The costs of carrying the alternate- 


* This formula, in the more general form of eq. (5), was first given by Rapp (Ref. 2). 

t The equation again is not strictly correct since the cost of the alternate route is not 
proportional to its offered load. However, we shall not use eq. (4) to evaluate the absolute 
cost, but only to determine its relative minimum with respect to n. For this purpose, the 
equation yields accurate results. 
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Fig. 5—Multihour trunk-group sizing for the network in Fig. 4. 


route loads in Hours 1 and 2 are drawn in the figure and are seen to in- 
tersect. To the left of the intersection point, the cost of the alternate route 
is determined by the Hour-1 load, and to the right, the cost is determined 
by the Hour-2 load. The maximization operator selects these portions 
of the curves; this is suggested by the shading in the diagram. The total 
cost is the sum of the alternate-route cost and the straight-line direct- 
route cost and is represented by the solid-line curve on the top of the 
diagram. The curve shows a minimum at the point n = no and has a 
discontinuous derivative at this point. The calculation of no for this 
example, then, differs radically from the calculation in our previous 
examples. Whereas before, the minimum-cost design was determined 
by requiring the load on the last trunk to be equal to a prescribed value, 
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here no is determined to be that value of n which equalizes the loads on 
the alternate route in the two hours. 

It is instructive to observe what would happen if a single-hour method 
were used in this example. If Hour 1 were chosen as the engineering base 
hour, the resulting high-usage group size would be nj, shown in the di- 
agram. Note that if n; trunks were installed, the alternate route would 
actually be busier in Hour 2 than in Hour 1, contrary to what was as- 
sumed. The total cost of the network after adding trunks to the alternate 
route to handle the Hour-2 load would be higher than at the optimum 
point no. Similarly, if Hour 2 were the selected engineering base hour, 
the resulting group size would be no. The alternate route would actually 
be busier in Hour 1, and again the total cost would be higher than opti- 
mum. 

Up to this point we have, for simplification, been treating the case 
where the alternate route has consisted simply of a single trunk group. 
Actually, of course, the alternate-route configuration is more compli- 
cated. Figure 6 shows one possible arrangement where the alternate route 
consists of a final group, a tandem switch, and a tandem-completing 
group. The background loads for the two hours are F'; and F2 for the 
final, S; and Sz for the tandem switch, and 7, and T's for the tandem- 
completing group. The cost per trunk of the final group is Cr and the 
cost per trunk of the tandem-completing group is Cr. We assume that 
the switching cost is proportional to the load and is equal to C's per CCS 
switched. The cost of each component of the alternate route again is 
determined by the maximum traffic offered to it. Thus, the engineering 
of the high-usage group requires three separate busy-hour comparisons. 
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Fig. 6—Typical alternate-route configuration. 
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The total cost is given by the formula 
Cr F, + a,B(n,a1) S, + a,B(n,a)) 
COST =—M 
7 os Fo + a2B(n,a2) Sot ea 
Cr Ti + a,B(n,a1) 
+—M 
y - Es + a2B(n,a2) 


| + Cs Max | 


+Cpn, (5) 


which is explained exactly like eq. (4). 

Figure 7 shows one possible example of the behavior of the cost curves 
as a function of n. In this example, the cost curves of the final group 
during Hours 1 and 2 intersect at a certain point. The switch cost curves 
also intersect, but at a different point. The tandem-completing cost is 
completely dominated by the Hour-1 load in this example. The sum of 
these three costs, plus the straight-line direct-route cost, yields the total 
cost curve shown at the bottom of the figure. In this example, the opti- 
mum design requires equalization of the switch loads in Hours 1 and 2. 
Of course, depending on the loads and costs, the minimum-cost point 
could instead have required equalizing of the final loads or of the tan- 
dem-completing loads. Alternatively, the minimum-cost network might 
not correspond to any of these “breakpoints” of the curve, but could lie 
on a smooth portion as in the single-hour case discussed previously.* 

It is important to observe that even in the case where the optimum 
solution does not correspond to a breakpoint, the engineering does not 
necessarily reduce to the single-hour ECCS method. To use the single- 
hour ECCS method, we first compute the “effective alternate-route cost 
per trunk” as C4 = Cr+ yCg + Cr. Then, using a single hour’s load, we 
determine n such that the load on the last trunk during that hour is equal 
to y divided by the cost ratio, C4/Cp. To see how the multihour method 
differs from this, let us assume that the final group is dominated by its 
Hour-i load, the switch is dominated by its Hour-j load, and the tan- 
dem-completing group is dominated by its Hour-k load. (We make this 
assumption to avoid having to worry about breakpoints in the cost curve.) 
Then, it follows from the differentiation of eq. (5) that the optimum value 
of n satisfies the equation 


La) + CsL (aj) + Lea) = Cp, (6) 


* Rather than attempting an exact minimization of eq. (5) in the manner we have de- 
scribed, Rapp adopted an approximate approach. He introduced the fictitious load, a, as 
a function of the parameters aj, a2, Ai, Ag, and then used it in the single-hour formula 
of eq. (3) to produce a trunk-group size.” 
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Fig. 7—Multihour trunk-group sizing for the network in Fig. 6. 


where a,, (m = 1,),k) is the load offered to the high-usage group during 
Hour m and L(a,,) is the load on the last trunk during that hour.* [This 
notation omits explicit indication that L(a,,) is a function of n.] Ifi = 
J =k, so that all components of the alternate route are busy at the same 
time, then L(a;) = L(a;) = L(a,) = L(a), and we can factor out this 
quantity. The equation becomes, in this case, 


* More precisely, L(am) = —(d/dn)|[a,B(n,am)]. 
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a6 


Ee Cea 


(7) 


where C4 = Cr + yCg + Cr. This is precisely the ECCS formula. [See 
eq. (3)]. 

Thus we see, as expected, that the single-hour ECCS method is opti- 
mum only when all components of the alternate route have the same 
busy-hour (or when the high-usage load is the same in all hours). If this 
is not the case, then it is invalid to lump the various components into one 
alternate-route cost and to choose only a single hour’s load for engi- 
neering the high-usage group. Equation (6) shows that, in general, the 
optimal sizing of the group involves the load on the last trunk in each 
of the three hours that are significant to the alternate route. 

An example is given in Fig. 8 which illustrates the difference between 
multihour and single-hour engineering in the case where the busy hours 
of the alternate-route legs do not depend upon the high-usage group size. 
In the example, a single high-usage group with offered loads of 300 CCS 
in the daytime and 200 CCS in the evening overflows to a final group with 
a known daytime busy hour and then to a tandem-completing group with 
a known evening busy hour. (The tandem is neglected for simplicity.) 
Since the busy hours of the alternate-route legs are fixed, eq. (6) gives 
the multihour solution for this example. The equation yields a trunk 
requirement of 11. Single-hour engineering produces a trunk require- 
ment of 13 for the daytime load and 9 for the evening load. It can be seen 
in the figure that the multihour network cost is significantly lower than 
both single-hour networks in this example. Note that daytime engi- 
neering over-sizes the high-usage group due to an overestimate of the 
tandem-completing cost, while evening engineering undersizes the group 
due to an underestimate of the final cost. 


IV. THE SIGNIFICANT-HOURS ALGORITHM 


An alternative procedure for engineering networks for more than one 
hour of traffic data, called the “significant- hours” method, has recently 
come into use in the Bell System. In this section, we describe this algo- 
rithm as applied to a two-level local network, and compare it to the 
multihour method discussed above. 

The significant-hours algorithm was devised to overcome shortcom- 
ings of the cluster-busy-hour ECCS approach caused by the fact that 
the various legs of the alternate-route path have busy hours different 
from that of the originating cluster. These shortcomings can be explained 
by considering the network in Fig. 9. In the cluster-busy-hour approach, 
the originating cluster-busy hour of each office is used for sizing all the 
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Fig. 8—Comparison of engineering methods for fixed alternate-route busy hours. 


originating groups in each office. Thus, group A-Z is sized for its A-office 
cluster-busy-hour load, group B-Z for its B-office cluster-busy-hour load, 
and group C-Z for its C-office cluster-busy-hour load. If offices A,B, and 
C are business-dominated offices, the groups A-Z, B-Z, and C-Z would 
be sized for their daytime business loads. If office Z is a residence- 
dominated office, however, the loads on these groups may peak in the 
evening. Since the groups are sized for their smaller daytime loads, they 
would overflow heavily in the evening, and all this overflow would be 
offered to the tandem-completing group T-Z. This effect has been ob- 
served in actual networks; in some cases, extremely high loads occur on 
certain tandem-completing groups, requiring great quantities of trunks 
and switching termination equipment. The problem is clearly caused 
by the exclusive attention to the originating portion of the alternate- 
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Fig. 9—Example of network. 


route path and the total neglect of the terminating portion. 

The significant-hours method solves this problem by giving equal 
treatment to both parts of the alternate-route path. For group A-Z, two 
significant hours are considered: the A-office originating cluster busy 
hour (defined above) and the Z-office terminating cluster busy hour 
(defined as that hour for which the total traffic terminating at office Z 
is maximum). Of these, the one for which the A-Z load is larger is chosen 
as the “control hour” for group A-Z. The group is then engineered for 
this load. By engineering the high-usage group for the larger of its sig- 
nificant loads, enough trunks are installed to eliminate the possibility 
of extremely heavy overflow in the busy hour of the tandem-completing 
leg. 

There still remain two problems with this method, however. The first 
is the fact that the actual busy hours of the final and tandem-completing 
groups of the engineered network are not necessarily the same as those 
that the significant-hour calculation predicts.* The second problem is 
that even if the significant hours are the right ones, in the sense that the 
alternate-route busy hours after engineering agree with the original as- 
sumptions, this method will over-engineer the group unless either (7) 
the significant hours are all the same, or (ii) the offered loads on the 
high-usage group are the same in these hours. 


* It should be emphasized that using observed final-group and tandem-completing- 
group busy hours, instead of the originating and terminating cluster busy hours, does not 
get around this problem. The observed busy hours depend on the previous network con- 
figuration. The point is that if a significant amount of noncoincidence of traffic loads exists, 
the busy hours of the newly engineered network will not agree with those assumed in the 
engineering, no matter how the hours are selected. 
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V. MULTIHOUR ENGINEERING OF A NETWORK 


In the discussion of the theory of multihour engineering surrounding 
eqs. (4) and (5), we considered only the sizing of a single high-usage trunk 
group in isolation. The background loads on the alternate route were 
given and assumed fixed. In a network, however, these background loads 
consist partly of overflows from other high-usage groups—groups which 
themselves have to be sized during the engineering process. For a net- 
work consisting of more than one group, therefore, the use of eq. (5) alone 
is insufficient, since it does not account for the interdependence between 
the high-usage groups that arises through their mutual effect on the 
background loads. 

The optimal sizing of a network consisting of N high-usage trunk 
groups in fact requires the minimization of a cost function of N dimen- 
sions, instead of the one-dimensional cost function of eq. (5). An analysis 
of this optimal approach has been carried out by W. B. Elsner of Bell 
Laboratories*. All of our initial network results, however, including those 
described in the remainder of the present paper, were obtained using 
a simple iterative approach. These initial results allowed us to demon- 
strate the feasibility of multihour engineering and to quantify the order 
of magnitude of the associated cost savings. These preliminary findings 
justified the effort by Elsner to develop the exact algorithm. 

We begin the iterative process by choosing initial sizes for every trunk 
group in the network.* This allows us to compute overflows from each 
high-usage group and thus to determine the total background loads 
which are offered to all alternate-route groups. We then size each high- 
usage group in turn by minimizing its one-dimensional cost function. 
The background loads used in each case consist of the first-routed loads 
plus the overflow from all other high-usage groups, that is, all high-usage 
groups except the overflow from the group being sized. After sizing every 
group once, the background loads that appear on the alternate-route 
groups differ from what they are at, the beginning and, hence, the engi- 
neering procedure is iterated, each pass consisting of the resizing of every 
high-usage group. This process continues until the iteration conver- 
ges. 

An essential aspect of this procedure is the fact that the background 
loads are updated immediately after the sizing of each group and before 
the sizing of the next group in sequence. The background loads play a 
very important part in the process of multihour engineering and it is 
necessary that the computed background loads be accurate if the proper 
sizing is to take place. If the updating is not done promptly, the back- 


* In obtaining the results that follow, we initialized each group to be numerically equal 
to the largest load on the group measured in erlangs. 
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ground loads used in sizing will differ from their true values, and 
misengineering of the network will result. If the updating is delayed to 
the end of each complete pass, for example (as was done in our first at- 
tempt), the iteration could even fail to converge. 

When the iteration does converge, the resulting network has the 
property that, if all other groups are held fixed, each individual group 
is sized to minimize cost. It is possible, however, that at convergence 
network cost could be further reduced by changing the sizes of two or 
more trunk groups simultaneously. This simple iterative approach also 
has the property that the solution to which it converges is not unique; 
depending upon the initial trunk values assumed, and the order in which 
the groups are sized, the solution network can vary. Both of these un- 
desirable properties of the iterative method are overcome with Elsner’s 
approach.? 


Vil. COMPUTER PROGRAM AND RESULTS 


A computer program incorporating the above iterative multihour 
procedure was written to design a network with the routing structure 
shown in Fig. 10. In this network, a single end office has a number of 
one-way outgoing high-usage trunk groups connected to other end of- 
fices. All high-usage groups overflow to a common final group and traffic 
reaches its terminating office via a one-way tandem-completing 
group. 

The program was run using the load data from three California end 
offices: Gardena, Compton, and Melrose. In each case, two hours of loads 
were employed. Hour 1 was a morning busy hour dominated by business 
traffic, and Hour 2 was an evening busy hour dominated by residential 
traffic. In the absence of actual trunk and switching cost data, the trunk 
cost of every group was assumed identical, equal to $1000 per trunk, and 
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Fig. 10—Network configuration for multihour engineering program. 
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the switching cost was assumed to be $62 per CCS (which yields a real- 
istic average switching-to-trunk-cost ratio). Both the switch and the 
tandem-completing groups were assumed to have zero background load. 
(This implies that the switch loads are identical to the final loads and 
that the tandem-completing busy hours are the same as their corre- 
sponding high-usage group’s busy hours.) 

Table I contains detailed results obtained for the Gardena office. The 
table shows the two hourly offered loads for each high-usage trunk group; 
it also shows the number of trunks and the resulting hourly overflows 
for the cluster-busy-hour network and the multihour network. 

From the totals at the bottom of the table, it can be seen that the total 
cluster load in Hour 1 is 6712 CCS while the total cluster load in Hour 
2 is 5154 CCS; clearly, Hour 1 is the cluster busy hour. On the other hand, 
for the network engineered for the cluster busy hour, the total overflow 
(i.e., the load offered to the final group) is 975 CCS in Hour 2 and only 
502 CCS in Hour 1. Here we see again the phenomenon of the final busy 
hour differing from the cluster busy hour. 

The reason for the high side-hour overflow can be seen by looking at 
the Hour-1 and Hour-2 overflow columns. The overflows in Hour 1 are 
quite uniform over all trunk groups; the trunk sizes are “matched” to 
the Hour-1 loads. In Hour 2, however, there is a great mismatch. A few 
trunk groups have very large overflows, the rest have virtually none. 
Three groups alone (14, 18, and 20), in fact, account for almost 60 percent 
of the total overflow in this hour. As can be seen, the pattern of overflow 
in the multihour network is more nearly balanced between the two hours; 
the total overflows in Hours 1 and 2 are nearly equal in this network. 

Table II compares the main characteristics of the busy-hour and 
multihour networks for Gardena. The numbers of final trunks shown 
are sufficient to guarantee a blocking probability of 0.01 for both hours. 
For the busy-hour network, for example, the final must be sized for the 
Hour-2 load since that load is larger. The total cost shown is the sum of 
the costs of the high-usage groups, final group, tandem switching, and 
tandem-completing groups. (The tandem-completing cost is an ap- 
proximation based upon the use of marginal capacity to determine the 
trunk requirement for each group.) 

With this simple model, the total cost of the multihour network is 
approximately 7 percent less than the cost of the single-hour network. 
Also significant is the fact that the switching cost is reduced 26 percent 
by using the multihour technique. Table III shows cost comparisons of 
all three of the networks studied. As can be seen, total network costs for 
the three cases decrease in the range of 5 to 11 percent and tandem- 
switching costs decrease in the range of 17 to 26 percent with the use of 
the multihour technique. (Trunk costs of $1000 and switching costs of 
$62/CCS were assumed in each case.) 
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Table | — Gardena network results 








Offered Load Cluster-Busy-Hour 
Trunk (CCS) Engineering Multihour Engineering 
Group Hour-1 Hour-2 Hour-1 Hour-2 
Overflow Overflow Overflow Overflow 


Hour! Hour2 Trunks (CCS) (CCS) Trunks (CCS) (CCS) 








1 60 140 3 10 62 4 4 42 
2 119 9 6 8 0 3 45 0 
3 82 20 4 10 0 4 10 0 
4 305 76 12 20 0 6 126 1 
5 30 0 2 5 0 0 30 0 
6 59 7 3 9 0 1 37 1 
7 102 56 5 10 1 4 19 3 
8 256 161 11 13 1 8 47 8 
9 366 230 15 15 0 12 46 4 
10 469 310 18 20 1 18 20 1 
11 115 115 5 15 15 5 15 15 
12 144 34 7 9 0 7 9 0 
13 206 335 9 13 80 10 7 61 
14 310 650 13 13 233 16 3 154 
15 284 319 12 13 25 12 14 25 
16 93 152 4 15 50 5 7 33 
17 17 24 1 5 10 1 5 10 
18 74 325 4 8 200 6 1 143 
19 102 158 5 10 37 5 10 37 
20 137 322 6 14 141 8 3 92 
21 222 247 9 18 28 10 11 18 
22 252 390 11 12 78 12 7 59 
23 445 194 17 21 0 17 21 0 
24 176 86 8 11 0 8 11 0 
25 83 29 4 ll 0 4 11 0 
26 98 21 5 9 0 5 9 0 
27 158 74 7 13 0 7 13 0 
28 124 36 6 10 0 6 10 0 
29 54 25 3 7 1 3 7 1 
30 38 1 2 8 0 2 8 0 
3l 31 17 2 5 1 2 5 1 
32 140 46 6 15 0 6 15 0 
33 96 30 5 8 0 5 8 0 
34 122 62 6 9 0 6 9 0 
30 163 57 7 15 0 th 15 0 
36 163 72 7 15 0 7 15 0 
37 296 238 12 17 5 12 17 5 
38 33 28 2 6 4 2 6 4 
39 240 3 10 16 0 10 16 0 
40 136 7 6 14 0 6 14 0 
41 54 4 3 7 0 3 7 0 
42 52 35 3 7 2 3 7 2 
43 206 9 9 13 0 9 13 0 
Totals 6712 5154 295 502 975 287 713 720 
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Table Il — Main characteristics of busy-hour and multihour 
networks for Gardena 

















Cluster 
Network Busy-Hour Multihour 
Characteristics Engineering Engineering 
High-usage trunks 295 287 
Final trunks 37 29 
Switching cost $60,450 $44,640 
Total cost $437,106 $405,315 
Table Ill — Cost comparisons for three networks 





Network Total Costs ($1000) Tandem Switching Costs ($1000) 
Office Busy-Hour Multihour % Decrease Busy-Hour Multihour % Decrease 





Gardena 437 405 | 60.4 44,7 26 
Compton 493 441 11 87.9 64.7 26 
Melrose 303 288 5 45.2 37.6 17 








Vil. CONTINUING WORK 


In this paper, we have described the basic theory of multihour engi- 
neering and have demonstrated, in a few simple networks, the potential 
savings it can bring about. A number of questions require answers for 
this technique to achieve acceptability for use in the field. These ques- 
tions include the following: the “hours” of data used for engineering may 
occur in different seasons of the year as well as different times of the day. 
How many hours should be included in the engineering of a network and 
how should these hours be determined? How is multihour engineering 
to be accomplished in a large-scale network with more than two levels 
in the hierarchy? How is trunk administration to be carried out in a 
multihour environment? What changes are required for two-way trunk 
groups? These questions and others have been the subject of intensive 
study at Bell Laboratories, and the answers will be reported on in future 
papers. 


Vill. CONCLUSIONS 


Multihour engineering is a technique that can provide significant 
benefits in the design of alternate-route traffic networks. In a comput- 
erized engineering environment, especially since automated data-col- 
lection methods make it possible to collect larger amounts of (and more 
accurate) traffic data, this technique should prove to be a realistic and 
preferable alternative to the older single-hour techniques. 
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Design of Quantizers for Real-Time Hadamard- 
Transform Coding of Pictures 


By F. W. MOUNTS, A. N. NETRAVALI, and B. PRASADA 
(Manuscript received June 29, 1976) 


A methodology is developed to obtain subjectively optimum quan- 
tizers for Hadamard-transformed still pictures. To exploit the per- 
ceptual redundancies that depend upon the local properties of the 
picture, a small block (2 X 2 X 2, horizontal-vertical-temporal) is used. 
A series of subjective tests was carried out to determine the visibility 
of impairment in the reconstructed picture when noise, which simu- 
lated the quantization noise, was added to the Hadamard coefficients 
in the transform domain (H-noise). A design procedure for quantizers 
was developed using these visibility functions. These quantizers min- 
imize the ‘mean-square subjective distortion” (MMSSD) due to quan- 
tization noise. The resulting picture quality and entropy were com- 
pared with that of Max-type quantizers which minimize the “mean- 
square error” (MMSE). This comparison indicates that the MMSSD 
quantizers based on subjective visibility of the quantization noise are 
less companded than the MMSE quantizers. Also for the same number 
of quantization levels, pictures coded with MMSSD quantizers have 
better quality and less entropy than the pictures coded with minimum 
mean-square quantizers. 


I. INTRODUCTION 


A methodology is developed in this paper to establish fidelity criteria 
that characterize human observers’ perception of noisy transform-coded 
pictures and to obtain optimum quantizers for the transform coefficients 
based on these fidelity criteria. The perceptual effects of impairments 
introduced in the transform domain are, in general, quite different from 
the impairments introduced in the picture domain. Our experiments, 
which are performed with the Hadamard transform of a stationary 
picture, determine the visibility of impairments in the reconstructed 
picture when noise (H-noise), which simulates the quantization noise, 
is added to a Hadamard coefficient. Functions that give the appropriate 
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subjective weighting of the quantization noise as a function of the 
quantity to be quantized are derived from these experiments. These 
functions, called visibility functions, are used in a systematic way to 
design quantizers for PCM and DPCM coding of the transform coeffi- 
cients. These quantizers are compared with minimum mean-square error 
(MMSE) quantizers, both in terms of picture quality and bit rates. 


1.1 Relationship to previous work 


Considerable attention has been paid to the transform domain in the 
recent work on picture coding.!-!4 Transform domain processing has 
several potential advantages. It produces less correlated (but not nec- 
essarily independent) coefficients. It redistributes the image energy so 
that a large amount of energy is packed in a few of the coefficients. 
Moreover, on inverse transformation at the receiver, both noise from 
quantization of coefficients and the channel errors get distributed over 
the block in a manner given by the inverse transform of a particular 
coefficient. 

A number of different transforms have been investigated; among them 
are: Karhunen-Loeve, Fourier, Hadamard, Haar, cosine, and slant 
transform. There have been several attempts®15.!6 to compare the various 
transforms and to find their relative merits for coding of pictures. Almost 
all of these comparisons have been with respect to the following three 
criteria: (i) the correlation between the coefficients, (ti) the mean-square 
approximation error caused by setting some of the coefficients to zero, 
and (iii) the computational complexity in obtaining transform coeffi- 
cients from picture elements (pels) and vice versa. Perceptual factors 
and the dependence of the picture quality on the particular transform 
and the block size have not received the attention they deserve. 

The irreversible processing of the transform coefficients, which de- 
termines the trade-off between picture quality and bit rate, has been 
performed in a number of ways; for example, (1) zonal sampling or 
masking, which drops some predetermined higher-order coefficients; 
(ii) threshold sampling, which drops those coefficients whose values are 
below a predetermined threshold (a certain amount of addressing in- 
formation must be sent in using this technique); (111) quantization of the 
coefficients—both amplitude (PCM) and differential amplitude (DPCM) 
quantization®!4 have been considered. Most of the work on quantization 
of the coefficients has centered around minimization of mean-square 
error as a criterion in designing the quantization characteristics. Several 
assumptions on the probability of the coefficients have been made, in- 
cluding the familiar gaussian case,!° to carry out this minimization. 
Exploitation of the psycho-visual properties of the viewer and the op- 
timization of the quantizers for the best subjective quality of the picture 
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has often been mentioned in the literature; however, no systematic 
methods are available for achieving these. 

Our work on obtaining the quantization characteristics may be com- 
pared with that reported by Landau and Slepian! and Tasto and Wintz.* 
For this reason, we give a brief review of their reports. In both only a 
single frame of picture data is used despite the fact that the quantization 
noise is more visible when a sequence of frames of the same scene is 
coded. 

Landau and Slepian considered both Karhunen-Loeve and Hadamard 
basis vectors for the linear transformation and found that the 
Karhunen-Loeve transformation required solution of an almost de- 
generate eigenvalue problem. They then used Hadamard transformation 
with a 4 X 4 block. The number of quantization levels given to each of 
the first ten coefficients was approximately proportional to the variance 
of that coefficient, and the last six coefficients (H;, to Hig) were dropped. 
The first coefficient was quantized by a 64-level uniform quantizer. 
Coefficients Hy through Ho were quantized with quantizers having a 
companding characteristic given by a function of the form y = kv‘x. 

Two arguments led Landau and Slepian to this quantization strategy. 
Firstly, since the variances of the lower coefficients are in general larger, 
coding them more accurately reduced the mean-square error. Secondly, 
the higher coefficients tend to be large in the busy regions of the pictures, 
where the viewers have more tolerance to amplitude errors. Thus, they 
used in an empirical way the consideration of a characteristic of the 
viewer as well as the statistics in the design of quantizers. They carried 
out over 100 experiments in which the decision levels and the repre- 
sentative levels of the quantizers were changed. However, since the 
number of choices is so large, their search could not be exhaustive and, 
therefore, their quantizers are the best only among those that they in- 
vestigated. 

Tasto and Wintz proposed an encoder using a 6 X 6 adaptive 
Karhunen-Loeve transform whose coefficients are quantized by what 
the authors call a “subjectively” optimized system of quantizers. This 
is done by first starting with a quantizer that minimizes the mean-square 
quantization error and then changing it by a trial-and-error procedure 
to obtain the “best” picture quality in the authors’ judgment. The “best” 
is again from among those encountered in the trial-and-error procedure. 
They also conducted subjective rating experiments to compare the 
performance of the minimum mean-square quantizers with the “best” 
quantizers. 


1.2 Basic objectives and approach 


Our basic objectives are to obtain fidelity criteria in the transform 
domain which incorporate psycho-visual properties, and to develop 
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systematic methods for the optimum design of coders based on these 
fidelity criteria. As mentioned above, perceptual properties of the human 
viewer have not been given sufficient importance in the transform-coding 
literature and, consequently, good models do not exist to explain the 
‘subjective effects of the quantization errors in the coefficients when the 
coefficients are inverse transformed to obtain the picture element 
(pel). 

In the pel domain, some efforts!”7-19 have been made to measure 
properties of human vision in psychophysical experiments and then 
utilize these to design coders. It is not easy to extend or utilize these 
techniques for the transform domain where we deal with blocks of pels 
instead of one pel at a time. Imperfect reproduction of coefficients of the 
block distributes distortion over the entire block upon inverse trans- 
formation. 

To take advantage of both the perceptual and statistical properties, 
some of the factors one has to study are: 


(t) Spread of the quantization error by inverse transformation. 
(it) Visibility of the quantization error in different coefficients. 
(tit) Statistical decorrelation. 
(iv) Probability distributions of the coefficients. 


In this paper, we do not attempt to solve this general problem but 
restrict ourselves to nonadaptive coding of stationary pictures using a 
2 X 2 X 2 (horizontal-vertical-temporal) Hadamard transform. Although 
a temporal structure of the block is not relevant for still pictures, it will 
be used in the next phase of our work which will treat coding of a se- 
quence of pictures. The Hadamard transformation has been chosen for 
its simplicity in implementation. The objective in choosing a small block 
is to exploit the perceptual redundancies which depend on local prop- 
erties of the picture. The small block ensures that the quantization noise 
can be placed in parts of the picture where it is least visible. However, 
it does result in some loss of coding efficiency on statistical grounds. To 
compensate at least partially for this, we also discuss the differential 
coding of the first transform coefficient H. 

In Fig. 1 the definition of the Hadamard coefficients for the block size 
2 X 2 X 2 is given. H is the sum of the element brightnesses within the 
block. H» is the sum of the line differences within the block. H3 is the 
planar difference. H4 is the sum of the element differences within the 
block. It may be noted that for stationary pictures, H5, Hg, H7, and Hg 
are all zero; further, any noise added to the first four coefficients gets 
repeated in the reconstructed signal at half the frame rate due to the 
block structure. As mentioned earlier, in coding frames of a single picture 
scene, the “nonmoving” noise patterns are, in general, less annoying than 
the moving noise patterns normally encountered in a television system 
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Fig. 1—Definition of Hadamard coefficients. The pel positions are A, B, C, D, E, F, G, 
H. The Hadamard coefficients are H1, Ho, H3, H4, Hs, He, H7, Hs. 


and, for this reason, a system was built to give a more realistic repre- 
sentation of television coding impairment. This system is described in 
Section II. 

Our method for determining the visibility functions of the noise in Ho, 
H3, and H, consists of the following. We add H-noise (which simulates 
the quantization noise) to a coefficient whenever its magnitude exceeds 
a threshold. This is done because each of these coefficients consists of 
difference quantities of pels and, therefore, may be expected to mask 
the noise as some function of their amplitude. For the DPCM coding of 
H, we add H-noise whenever the magnitude of the difference of H,; from 
its previous block value is higher than a threshold. Again, this difference 
of H, can be taken as a measure of signal busyness. The effect of this 
H-noise impairment on the picture is then compared by the subject in 
an A-B test with simple additive white noise impairment of the picture. 
This method of judging pictures is similar to the one used by Candy and 
Bosworth.2° The experimental method is discussed in detail in Section 
II. 


71.3 Summary of results 


The visibility functions for the following conditions have been mea- 
sured: H»-noise as a function of |H2|; H3-noise as a function |H3|; H4- 
noise as a function of |H4|; and H-noise as a function of | AH,|, where 
AH, is the adjacent block difference in the horizontal direction. 

The study of these visibility functions indicates that Hs is the least 
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important coefficient and can be dropped entirely with little impairment. 
PCM quantizers with minimum mean-square subjective distortion 
(MMSSD) have been designed for Hz and H, coefficients, and an MMSSD 
DPCM quantizer has been designed for H, using the corresponding vis- 
ibility functions as the fidelity criteria. These quantizers have been 
implemented and have been compared in subjective tests with the cor- 
responding quantizers optimized with respect to the minimum mean- 
square error criteria. Details of this approach and the results are given 
in the subsequent sections. 


ll. EXPERIMENTAL SYSTEM 


The experimental system described in this section has been designed 
with considerable flexibility as a vehicle for future research. The system 
has real-time capabilities for adaptive and nonadaptive Hadamard 
transform coding of a 2 X 2 X 2 block of pels. 


2.1 System block diagram 


A block diagram of the experimental system is shown in Fig. 2. The 
video signal is generated by a vidicon camera scanned with 271 lines 
interlaced 2:1. The video signal has a bandwidth of 1 MHz and is sampled 
at the Nyquist rate. Each picture sample is PCM encoded with amplitude 
accuracy of 8 bits per pel. 

A frame memory is incorporated in the system to accommodate the 
transform block. Alternate frames of the digitized pictures, say the odd 
frames, are stored in the frame memory via data select switch 1. Memory 
1 consists of two line delays and four small delays for linking the data 
from the present and previous frames. It ensures, during even frames, 
simultaneous presentation of all the elements from the two frames that 
comprise the data block to the Hadamard transform logic. It may be 
noted that the system is designed for spatially overlapped block pro- 
cessing. The output corresponding to nonoverlapping blocks is selected 
by memory 2 for decoding and experimentation. The spatially overlap- 
ping blocks are suitable for the study of various kinds of predictive en- 
coders. This facility is also very useful for a flicker-free display of the 
Hadamard coefficients on a television screen. 

The Hadamard transform circuit is a serial and parallel combination 
of adders and subtracters to implement the canonic forms shown in Fig. 
1. In the processor circuit, the magnitudes of all the coefficients are 
rounded off to the eight most significant bits, which are used for further 
processing. This rounding off does not produce any visible impairment 
on inverse transformation. Capabilities exist in the processor circuit to 
insert eight independent quantizers, one for each coefficient. 

For the subjective experiments, the processor circuit permits two 
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Fig. 2—System block diagram. 


DIGITAL-TO- 
ANALOG 
CONVERTER 
















LP 
FILTER 










SUBJECT’S 
CONTROL 
EQUIPMENT 


RANDOM 
NOISE 
GENERATOR 






modes of operation controlled by the A/B switch. The details of the 
subsystem are shown in Fig. 3. The A/B switch is under the control of 
the subject. In the A mode, unimpaired coefficients are fed to the in- 
verse-transform circuit. This provides the original picture in the re- 
constructed signal domain. In the B mode, a controlled amount of 
pseudo-random noise is added to one of the coefficients only when the 
magnitude of a control signal (which is the coefficient itself in this di- 
agram) exceeds some reference threshold. This noise, which we call the 
H-noise, is generated at the sample rate by an 8-bit pseudo-random 
generator having a period of 2!5 words and which is not synchronized 
with the line or the frame rate of the picture. An amplitude limiter 
controls the magnitude of the noise to the level set by the experimenter. 
The sign bit for the noise word is obtained from the output of a white 
noise source, and has equal probability of being a “0” or a “1”. 

Since the addition of pseudo-random noise results in doubling of the 
maximum amplitude of the noisy coefficient, the sum of the coefficient 
and noise, and the other coefficients, are divided by 2 prior to inverse 
transformation to prevent overload. 

The inverse transformation network is similar to the transformation 
network and is used to reconstruct simultaneously all of the pels of the 
block. 

It may be recalled that the alternate frames (odd frames) of the input 
are stored in the frame memory via data select switch 1. The recon- 
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Fig. 3—Noise adding circuits. Coefficients are divided by 2 to prevent overload in the 
inverse transform function after noise addition. 
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Fig. 4—Original picture used for subjective tests. 


structed pels corresponding to the even frames are stored in the proper 
time slots in the frame memory by data select switch 1. Thus, the frame 
memory contains both processed and unprocessed data and is utilized 
fully. Data select switch 2 ensures that the reconstructed pels corre- 
sponding to the even frames that are stored in the frame memory are fed 
to the digital-to-analog converter in the proper time sequence. 

The original picture used for the subjective tests is shown in Fig. 4. 
The scanned and filtered version (by a 1-MHz Picturephone® filter) is 
shown in Fig. 5a. Figure 6 shows the picture of the coefficients using 
overlapping blocks. Figure 6a shows coefficient H,, which is essentially 
a “block-low-pass-filtered” version of the picture and preserves much 
of the picture information. On the other hand, Figs. 6b (H2 coefficient), 
6c (H3 coefficient), and 6d (H4 coefficient) show a variety of edge in- 
formation. 


2.2 Experimental details 


The experimental setup for determining a visibility function is shown 
in the simplified block diagram in Fig. 7. The experimenter adds H-noise 
to a selected coefficient whenever the absolute value of the coefficient 
exceeds a threshold. The amount of noise and the threshold are varied. 
This is presented as condition B to the subject. Condition A is the un- 
impaired picture plus white noise. By turning an attenuator knob, the 
subject can control the amount of white noise added to the unimpaired 
picture. He can switch between conditions A and B by the A/B switch 
provided. An experiment consists of the subject changing the attenuator 
until he finds the pictures in the switch positions A and B to be subjec- 
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(c) (d) 


Fig. 5—Filtered test picture and H-noise added pictures. (a) Filtered test picture (1-MHz 
Picturephone® filter). (b) Picture with noise added to H;. (c) Picture with noise added 
to H». (d) Picture with noise added to H4. 


tively equivalent. The subject can switch between A and B conditions 
as often as he likes and can look at the test conditions as long as he likes. 
When he arrives at the subjective equivalence, he gives the attenuator 
reading to the experimenters on an intercom. He is then given the next 
test condition. In one sitting, a subject makes 28 judgments of which the 
first four are considered as training. The remaining 24 are recorded as 
data. The experiment is also characterized by the following: 


(t) The picture has 271 lines, interlaced 2:1 at 30 frames per second. 
(ti) The visible portion of the picture is about 138 cm X 12 cm. 
(iii) High light brightness is 74 foot-lamberts. 
(iv) Low light brightness is 4 foot-lamberts. 

(v) Room illumination is 57 foot-candles. 


The scan lines of the Conrac monitor were broadened to correspond to 
the Picturephone display tube. Subjects were seated at a distance of 
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(c) (d) 


Fig. 6—Pictures of coefficients. (a) Picture of H, coefficient (with no output filter). (b) 
Picture of He coefficient (with output filter). (c) Picture of H3 coefficient (with output 
filter). (d) Picture of H4 coefficient (with output filter). 


about 80 cm from the monitor. All of the six subjects used had experience 
in judging coded television pictures. 


lil. TEST DATA AND ANALYSIS 


Results of a typical subjective test are shown in Fig. 8. In this case, the 
absolute value of H4 was compared to a threshold and the noise (H-noise) 
was added to H. In this figure, H4-noise is plotted in dB on the X-axis 
and the “equivalent white noise” is plotted on the Y-axis. Each datapoint 
is an average of the readings obtained from six subjects. Under the as- 
sumption that the equivalent white noise (Vy) is proportional to the 
H-noise (V77), the results for each threshold should fall on a 45-degree 
straight line. The lines drawn in Fig. 8 are the best unity-slope straight 
lines obtained by the least square fitting to the datapoints. Figures 9 and 
10 show similar data for Hz, and Hy, respectively. In the case of H,, | AH;| 
is compared to a threshold. Notice that in each case the quantity that 
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Fig. 8—Plot of “equivalent white noise” vs H4 noise for different thresholds on |H4]. 


is compared to a threshold is a measure of busyness of the picture in a 
local area. The locations in the picture where noise is added and its ap- 
pearance are dependent upon the quantity that is compared to a 
threshold and the coefficient to which the noise is added. 

The pictures with H-noise impairments are shown in Fig. 5. The H- 
noise added in each of the three pictures has a peak value of 100 units* 
(signal range is 0 to 255 units). In Fig. 5c, noise is added to H¢ in all blocks 
which have |H»| more than five units, whereas in Fig. 5d, noise is added 
to all blocks in which |H4| is more than five units. While H is the line 
difference, Fig. 5c has noise whenever an edge has a sufficiently large 
component along the horizontal direction, whereas Fig. 5d has noise 
whenever an edge has a sufficiently large component along the vertical 
direction. Also notice the difference in the appearance of the noise. 
Ho»-noise is much more noticeable than H4-noise. Figure 5b shows noise 


* This is much more than the noise used for any test condition, but has been used to 
demonstrate the effects in a photograph. 
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Fig. 9—Plot of “equivalent white noise” vs H2 noise for different thresholds on |H9|. 
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of 100 units added to H; whenever |AHj,| is more than 10 units. Here 
again, noise gets added to all blocks having horizontal interblock edges. 
Also the pattern generated by H, noise is much more objectionable than 
the pattern generated either by H2-noise or H4-noise. 

The relationship of proportionality between the equivalent white noise 
and the H-noise for a threshold z is written in the form: 


Vw = F(z) Vu, (1) 


where z can take on the value of |Ho|, |H3|, |H4|, or |AH1|. The constant 
of proportionality F(z) is the equivalent white noise power when a unit 
H-noise is added to the particular coefficient for all blocks of the picture, 
where the magnitude of the corresponding coefficient (| AH,| in the case 
of H,) is greater than or equal to the threshold z. We next assume the 
additivity of the equivalent white noise power with respect to the coef- 
ficient value; i.e., if the equivalent white noise power when a unit of 
H-noise is added to Hz and T, < |H2| < T2is V1, and the equivalent 
white noise power when a unit of H-noise is added to Hz and T»2 S |H2| 
<73(T; < T2< T3) is V2, then the equivalent white noise power when 
a unit of H-noise is added to Ha and T; S |He| < T3 is (Vu1 + Vuz2). 
Under this assumption, F(z) can be written as an elemental sum of the 
equivalent white noise powers. Thus, 


F(z) = f * f(x)dx, (2) 


where f(z) is called the visibility function. 

Using this procedure, visibility functions were computed. They are 
shown in Fig. 11. Notice that we have assumed that the occurrences of 
positive and negative coefficients (AH,, He, H3, H4) are similar, and the 
noise visibility does not depend upon the sign of the coefficient. This 
results in the visibility functions being symmetrical about zero. The value 
of the visibility function shows the relative importance of the various 
transform coefficients. The larger the value, the more important is the 
coefficient. In general, the visibility functions decrease as a function of 
their arguments. This is a combined effect of several factors, such as (2) 
the decrease in the number of blocks having large coefficient values 
(| AH;,| in the case of H,), (ii) the dependence of the perception of noise 
on the magnitudes of the coefficients (which correspond to the sharpness 
of the boundary in the pel domain), and (ii) the contextual importance 
of the specific regions of the picture. 

Psycho-visual techniques which measure the detectability of per- 
turbations in the neighborhood of edges,!®:!-23 and the just noticeable 
differences in the amplitudes of edges have been widely applied to DPCM 
coding. Since these deal with over-simplified stimuli and surround and 
are almost always detection experiments, their use in picture coding may 
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not always result in better coders. In any case, these techniques cannot 
be easily applied in the transform domain because we are not dealing 
with the single pels but with blocks comprising pels from more than a 
single line and frame. Also, the perturbations must be introduced in the 
transform coefficient, whereas the annoyance to the perturbations must 
be judged in the pel domain. 

Our approach, which obtains the visibility functions as outlined above, 
has the following limitations: 

(i) Since the visibility functions are tied to the picture content, they 
admittedly vary from picture to picture, especially if the picture content 
is changed significantly. They also depend upon the class of viewers and 
the viewing conditions. Thus, any optimization based on the visibility 
functions is strictly applicable to a restricted situation. This is under- 
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Fig. 11—Plot of visibility functions. Notation H;, H; indicates noise added to H;, when 
H; is thresholded. The H; threshold level is shown on the X-axis. 
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standable since the human perception does indeed vary with the picture 
content, the viewing conditions, and the particular viewer. We demon- 
strate in Section IV that the results we obtained using these visibility 
functions are not overly sensitive to the picture content and are rea- 
sonable for a class of pictures rather than a particular picture. 

(it) The simulation of the quantization noise by the H-noise is fairly 
accurate for Ho, Hs, and H4. However, in the case of H;, DPCM tech- 
niques are used to code AH. The changes in the appearance of noise as 
a function of threshold are not completely reflected in the measurements; 
i.e., while the noise that is added to H does look like granular noise at 
low thresholds, it does not look like slope overload at high thresholds. 
Also, at high thresholds, the noise is added to fewer blocks in the picture 
and the appearance of such an impaired picture is different from the 
appearance of a white noise impaired picture. Therefore, in some cases 
subjective equivalence is hard to achieve. 

(iii) It would be better if the perceptual, statistical, and contextual 
effects were explicit in the visibility function and could be controlled 
separately. Unfortunately, such is not the case. 

(iv) It is seen from eq. (2) that the process of obtaining the visibility 
function involves differentiation of the data, which is known to introduce 
some noise. By adding H-noise to a coefficient when the quantity to be 
compared to a threshold is within a small range of values, it is possible 
to avoid this differentiation. 


IV. RESULTS 


In this section, we present certain conclusions drawn from the visibility 
functions and then describe their application to the design of quantizers 
for the coefficients. Visibility functions shown in Fig. 11 clearly show 
the relative importance of various coefficients. H, is the most important, 
Ho is the next, followed by H4, and H3 is the least important. The visi- 
bility of H noise depends upon the patterns associated with a particular 
coefficient. These patterns depend upon the inverse transform and are 
shown in Fig. 12 for Ho, H3, and H4, respectively. In each case, noise of 
a given amount is added to one of the coefficients, and the background 
is assumed to be flat. The higher the spatial frequency of the pattern, 
the lower the visibility of the noise. Thus, H2 noise is more visible than 
H 4 noise because the interlace gives the H2 noise pattern lower spatial 
frequency than the H4 noise pattern. 


4.1 Visibility of H, noise 


An experiment was performed to utilize the well-known property of 
the human eye that the brightness discrimination decreases as the 
brightness level increases, called Weber’s law in the psycho-visual lit- 


QUANTIZER DESIGN 37 





Fig. 12—Pictures of noise patterns for coefficients. (a) Picture of noise added to H on 
a flat background. (b) Picture of noise added to H2 on a flat background. (c) Picture of 
noise aoded to Hs on a flat background. (d) Picture of noise added to H, ona flat back- 
ground. 


erature?4-°6 (see Ref. 27 for a recent application of Weber’s law for pic- 
ture coding in the pel domain). 

In this experiment, noise was added to H as a function of H, since 
H corresponds to the average brightness in the block. The results of this 
subjective experiment showed large variations from observer to observer. 
When the data for the observers was averaged, there was no significant 
variation in the visibility of noise as a function of H;. This could be due 
to the following: (i) If the gamma of the display tube used was not unity, 
it would have partially compensated for Weber’s law effects. (ii) We were 
working with the head and shoulders view of a person. In general, for such 
a picture, the highlights are on the forehead or the cheek of the person. 
These regions are contextually very important causing the visibility of 
noise to be high. (ii) The picture we used was such that the low-light 
areas had more spatial detail than the highlight areas; thus, the latter 
two effects may have compensated for the Weber’s law. 
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Measurements of the gamma of the monitor indicate that the visibility 
function in this case cannot be fully explained on the basis of the com- 
pensation of Weber’s law by the gamma of the display tube. It seems that 
at least for this class of pictures, namely the head and shoulders view of 
a person, the advantage that could be gained by the utilization of Weber’s 
law is compensated for by the other effects. 


4.2 “Frozen” vs “unfrozen” noise visibility 


It may be recalled that the experiment on visibility was done using 
a block size of 2 X 2 X 2. In this case, any noise added to the first four 
coefficients remained unchanged for two frames. The noise in the coef- 
ficients in this case may be called the “‘frozen” noise because it remains 
unchanged for two frame periods. An experiment was performed to de- 
termine the visibility functions for H4 coefficient for a block size of 2 x 
2 (horizontal-vertical). Since all the experiments have been carried out 
with a stationary picture, the only difference between the experiment 
with the block size of 2 X 2 and a block size of 2 X 2 X 2 is the coefficient 
noise. For the block size of 2 X 2, the coefficient noise changes from frame 
to frame and is called ‘‘unfrozen” noise. Figure 13 shows the visibility 
functions for H4 with “frozen” and “unfrozen” noises. Although the 
visibility functions of “unfrozen” noise are generally a little lower than 
that of “frozen” noise, due to lower temporal frequency, the differences 
are small. 

In Section III, it was mentioned that the visibility functions can be 
used as fidelity criteria for the design of quantizers. We describe below 
how these results are used to design quantizers. 


4.3 pcm coding of H2, H3, and H, 


It is assumed that little interaction exists between the Hadamard 
coefficients, so that the quantization transfer characteristics for the 
coefficients can be obtained independent of each other. It is recalled from 
Fig. 11 that H3 was the least important coefficient and, therefore, it was 
decided to drop the transmission of H3 altogether. 

Minimum mean-square error quantizers are obtained by minimizing 
the mean-squared quantization error. If N is the number of levels, and 
Pu,(-) is the histogram for | H;,|, then we minimize the distortion D given 
by 


N Xj+1 
De (|H,| — Y;)? Px, (|Hx|) d(| Ae), k = 2,3,4, (8) 
J=1 j 
with respect to {Xj}, 7 = 2,..., Nand {Y;},j =1,..., N. This gives us the 
well-known Max quantizer.22 MMSE quantizers are obtained by 
weighting the quantization error according to the frequency of its oc- 
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Fig. 18— Visibility function for “frozen” and “unfrozen” H,4 noise. 


currence. Minimum mean-square subjective distortion quantizers, on 
the other hand, weight the quantization error according to its subjective 
visibility. This can be achieved by substituting f7,(-) for P7,,(-) in the 
expression for the distortion. The term fy, (-) is the visibility function 
for the coefficient H;. Standard programming techniques were used to 
minimize the distortion D in both cases. 

The histograms for | AH,|, | H»], |H3|, and |H4| are shown in Fig. 14. 
In general, these decrease faster than the visibility functions. This is 
exemplified in Fig. 15 in which the histogram and the visibility function 
for AH, are plotted with the same scale on the X -axis. We shall see later 
that this fact results in larger companding of the MMSE quantizers than 
the MMSSD quantizers and, consequently, poor reproduction of busy 
areas of a picture. 

Typical quantizer characteristics are shown in Fig. 16. The MMSE 
quantizer is more companded than the MMSD quantizer. Note also that 
the dynamic range of the MMSE quantizer is smaller. The performance 
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Fig. 14—Histograms of coefficients. All the histograms are assumed to be symmetric 
about zero. 


of these two types of quantizers (MMSSD and MMSE) was compared in 
an A-B test with different numbers of levels. Figure 17 shows the results 
of such a test for the coefficient 14. In this test, MMSE quantizers with 
levels 5 to 8 were compared in terms of picture quality with MMSSD 
quantizers with levels 3 to 9 using a random pairing by six skilled 
subjects. The numbers in the table indicate the percentage of observers 
who preferred the MMSSD quantizers over the MMSE quantizers. The 
picture coded with the 5-level MMSSD quantizers was preferred by 100 
percent of the subjects over the 6-level MMSE quantizer. Figure 18 shows 
similar comparisons for the quantization of Hy. Here again, for the same 
number of levels, the picture quality using the MMSSD quantizers is al- 
ways better than using MMSE quantizers. Moreover, picture quality using 
the 6-level MMSSD quantizer and 7-level MMSE quantizer is equivalent. 
Figures 19 and 20 show the entropy of the quantized output using both 
the MMSSD and the MMSE quantizers having levels 3 to 8 for H4 and Ho, 
respectively. In the case of H4, the difference between the entropies of 
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Fig. 15—Comparison of probability density (P) and visibility (V) for H noise. 


the output of the MMSE and the MMSSD quantizer for the same number 
of levels is about 0.2 bit. Since the picture quality with the 7-level visi- 
bility quantizer was better than with the 8-level MMSE quantizer, the 
gain by the use of the MMSSD quantizer is of the order of 0.5 to 0.6 bit 
for the transmission of H4. Similar remarks can be made about the 
quantization of Ho. 

MMSSD QUANTIZER (8 LEVELS) 
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Fig. 16—Typical (MMSSD and MMSE) quantizer characteristics for H4 coefficient. 
Notation ie implies that all input levels between x and y (including x,y) are represented 
as Zz. x y 
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Fig. 17—Comparison of picture quality of MMSSD and MMSE quantizers of different 
levels for H4. 


4.4 Coding of H, 


Unlike the Ho, H3, and H4 coefficients, the H, coefficient is not a 
difference signal. It represents the average brightness within the block 
and thus carries the low-frequency information which should be coded 
relatively precisely. Uniform PCM coding of H, requires 7 to 8 bits for 
good picture quality. As mentioned before, efforts to compand the PCM 
quantizer by using the Weber’s law effect were not very successful. 
Therefore, it was decided to DPCM encode H). Since the block size used 
is small, there is substantial correlation between the H values of adja- 
cent blocks. This was exploited by using a DPCM coding of H, with ho- 
rizontally adjacent blocks for prediction. The quantizers for such a DPCM 
coder are obtained from the visibility function of H; under the control 
of |AH,| in a manner similar to the above by minimizing the mean- 
square subjective distortion due to the quantization noise. The resulting 
quantizer scales are companded due to the monotonic decrease of the 
visibility function with respect to AH), as shown in Fig. 11. Quantizer 
scales have also been obtained by minimizing the mean-square quanti- 
zation error.* As noted before, MMSE quantizer scales are more com- 
panded and have less dynamic range compared to the MMSSD quantized 
scales. Using these two types of scales, experiments have been performed 


* ~ * Although h the visibility function and the histograms are obtained from the difference 
signal | AH if and the quantity that is quantized is the differential signal (i.e., the difference 
between the present H, and the coded value of H) from the previous block), it is expected 
that the quantizer characteristics will not change appreciably by using difference instead 
of the differential signal in eq. (3). 
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Fig. 18—-Comparison of picture quality of MMSSD and MMSE quantizers of different 
levels for Ho. 
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Fig. 19—Plots of entropy of outputs of MMSSD and MMSE quantizers of different levels 
for Hyg. 
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f Hie. 20—Plots of entropy of outputs of MMSSD and MMSE quantizers of different levels 
or i119. 

to compare the picture quality for the same number of levels. The results 
of such a comparison are shown in Fig. 21. It is seen that, for the same 
number of levels, all the subjects preferred the picture coded with the 
MMSSD quantizers over the picture coded with MMSE quantizers. 
Moreover, picture quality using a 24-level MMSSD quantizer is equivalent 
to the picture quality using a 30-level MMSE quantizer. Entropies of the 
quantized signal with MMSSD and MMSE quantizers of different levels 
are shown in Fig. 22. Here again, visibility quantizers perform better than 
MMSE quantizers by about 0.35 bit per block for the same number of 
levels. Picture quality using a 24-level MMSSD quantizer can be produced 
by an MMSE quantizer with an increase in entropy of 0.6 bit per block. 
It is worth noting tht, due to the DPCM coding of H, the bits required 
for H, could be almost halved. However, 7, still remains more important 
than H» and H, and requires more bits for satisfactory transmission. 
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Fig. 21—Comparison of picture quality of MMSSD and MMSE quantizers of different 
levels for DPCM coding of H. 


4.5 Combined quantization of all coefficients 


Combined quantization of all the coefficients requires investigation 
of the optimal number of quantizer levels to be given to each one of them. 
In the case of gaussian random vectors using Karhunen-Loeve trans- 
formation and mean-square error criterion, optimal bit allocation for 
the various transform coefficients is well known.29:29 However, in our 
case, none of these assumptions are strictly valid. In fact, our assumption 
that the optimum quantizer characteristics for different coefficients can 
be obtained independently is not strictly true and, for this reason, we 
tried to evaluate the picture quality by quantizing all the coefficients. 
By trial and error, a near-perfect picture was produced by using 36, 13, 
and 7 quantization levels for AH,, He, and H, respectively, and by 
dropping H3. This resulted in a total entropy of about 2.17 bits per pel. 
In single-frame photographic reproduction, no difference could be ob- 
served between the coded picture and the low-resolution original shown 
in Fig. 5a. Several other “head and shoulders” type of pictures were coded 
using the same combination of levels. Although, in each case the picture 
appeared to have a reasonable quality, the visibility of the quantization 
and the resulting picture quality varied slightly. This implies that the 
quantizers we obtained by optimizing the visibility of the quantization 
noise for one particular picture were not overly sensitive to variation in 
picture content. 


V. CONCLUSIONS 


A systematic method for quantizing Hadamard coefficients has been 
given. This method gives the best quantizers in a subjective and proba- 
bilistic sense. We have compared the resulting quantizers with MMSE 
quantizers and found the MMSSD quantizers to be better both in terms 
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Fig. 22—Plots of entropy of outputs of MMSSD and MMSE quantizers of different levels 
for DPCM coding of H. 


of the subjective picture quality and entropy. We do not imply that there 
are no better quantizer than the MMSSD quantizer, since by taking many 
other factors into consideration, one could come up with a better 
quantizer. We do find that the minimum visibility quantizers are opti- 
mum with respect to our model and the approach used for weighting the 
quantization noise. 

Investigations are in progress for adaptive and predictive coding of 
the coefficients; our findings will be reported in a future paper.?! 
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Design criteria to minimize modal dispersion have been found for 
a broad class of practical, multimode, circular-symmetric, isotropic, 
optical fibers having any numerical aperture and any profile dispersion 
(which is a function of the derivative of the index with respect to the 
wavelength). The impulse-response width of these fibers, the rms width 
of the impulse response, the optimum profiles to minimize those widths, 
and the sensitivity to profile departures from ideal are found to be 
surprisingly simple closed-form generalizations of previous results that 
are mostly applicable to fibers with small numerical aperture and 
constant profile dispersion. The minimum impulse-response width of 
the optimized fiber is a function only of its numerical aperture and 
consequently is independent of the index profile and of the profile 
dispersion. 


|. INTRODUCTION 


Circular-symmetric, multimode, optical fibers intended for large 
communication capacity must have low modal dispersion and this is 
achievable by the quasi-complete equalization of the group velocities 
of all modes! (or rays). This equalization depends critically both on the 
refractive-index profile and on the profile dispersion of the fiber. The 
profile dispersion is defined in Section II, but here it is enough to know 
that it.is related to the derivative of the index with respect to the wave- 
length. . | . »t % 7 

To understand better the objectives of this paper, let us first review 
some recent evolution of thoughts linking the index profile and: the 
profile dispersion of a fiber to the pulse broadening caused by modal 
dispersion. 
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Gloge and Marcatili? showed that if the numerical aperture (NA) of 
the fiber is arbitrary but the profile dispersion is negligible, there is a 
family of fibers—for which the dielectric constant profiles decrease ra- 
dially according to power laws—that is important for two reasons. The 
first reason is that the family encompasses a wide variety of easy-to-make 
fibers (step-index, quasi-parabolic, etc.) possessing the unique property 
that the group velocity of each mode is a function only of its propagation 
constant; this drastically simplifies the analysis. The second and more 
important point is that for an almost parabolic power law of the dielectric 
profile, a fiber with small NA has the very narrow impulse response 
needed for high-speed communication. 

Olshansky and Keck* extended these results in a very important way 
by showing that if the profile dispersion is constant across the core, 
narrow impulse response is achievable in small NA fibers by a simple 
modification of the exponent of the dielectric-constant profile’s power 
law. 

In many cases, though, the two requirements—smallness of NA and 
constancy of profile dispersion—are not satisfied. For example, to in- 
crease the coupling efficiency to incoherent sources and to decrease 
microbending losses,’ fibers with large NAs® are being made. They are 
heavily doped and, particularly if the doping element is boron, the profile 
dispersion may not be a constant®’ as a function of the radius. Similar 
lack of constancy may occur in fibers that are doped with several mate- 
rials for the purpose of improving optical or mechanical properties.® 
Arnaud and Flemming?:!° have calculated the impulse response for these 
fibers, treating the variable portions of the profile dispersion as a per- 
turbation. Using a numerical method Arnaud!! has also calculated the 
pulse spreading in a multimode planar fiber with arbitrary index profile 
and profile dispersion. 

In this paper, we extend the previous results by finding, within the 
WKB approximation, a surprisingly simple closed-form description of 
the modal dispersion in a broad class of circular-symmetric isotropic 
fibers which have arbitrarily large NA and arbitrary profile disper- 
sion. 

The gist of our paper is in Sections IJ and III. In Section II, the profiles 
of the fibers belonging to the group are defined and their impulse-re- 
sponse widths are calculated. In Section III, the optimum index profile 
required to minimize the impulse-response width is determined together 
with the sensitivity of this response to departures of the index from 
optimum. The rms of the impulse response is the subject of Section IV. 
In Section V, some approximate results about the influence of index 
profiles on the minimization of the rms width of the impulse response 
are derived and conclusions are drawn in Section VI. 
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ll. FIBER PROFILE AND WIDTH OF ITS IMPULSE RESPONSE 


We start by looking for the dielectric-constant profile of a circular 
symmetric isotropic fiber such that, as in Ref. 2, the group velocity of 
each mode is only a function of its propagation constant. In the process, 
we will also find the width of the impulse response of that fiber. 

The initial point is a WKB approximation? that relates the propagation 
constant (A) of a mode characterized by the radial and azimuthal wave 
numbers u and » to the free-space propagation constant k = 27/X, the 
refractive index n(r, ) of the fiber and the radial coordinate r, via the 


integral 
1 ro dr 
b= f pv, (1) 
Tv ry r 
where 
p(r, \) = V (k2n2 — B2)r2 — v2 (2) 


and r, and ry are two neighboring turning points that make the radical 
zero and between which most of the field of the mode is concentrated. 
It is useful to redefine 


n2=ni(1 —F) (3) 

6? = k?n7(1 — B), (4) 

where n, is the index on axis and the profile function F(r, A) is zero on 
axis, is an arbitrary function of r and \ within the core (r S a), and is 
2A(A), a function only of ) in the cladding (r = a). Similarly, the mode 
parameter B varies between zero for the lowest-order mode and 2A()) 


for the modes whose phase velocities coincide with that of a plane wave 
in the cladding.* With these definitions, the radical (2) becomes 


p= V(knyr)?(B — F) — v?. (5) 


The group velocity of a mode (or ray) is introduced by taking the de- 
rivative of both sides of (1) with respect to the free-space wavelength 


r, 
re n, (dB Dp iE 

Bie Se ae i) ‘dr =0, 6 
S| ( aeTe, ( 2) 1p. (6) 

where 

A dn, 

N,= t= 7 
: na ( ny a) 7) 


* Similar but not identical profile function and mode parameter have been introduced 
previously in the literature.!2 
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is the group index on axis and 


ny \ OF 
p(r, A) N, Fan (8) 
is a generalized version of the profile dispersion parameter introduced 
in Ref. 12. 
The derivative dB/d) in (6) can be expressed in terms of the group 
delay t of the mode by taking square roots and derivatives on both sides 
of (4). The result is 


Cl ane rv dB 
=—=V1-B+— 
dkn, T a VI-Bdy @) 





in which T, the flight time on axis, that is, the delay of a plane wave in 
a medium of group index N, and length L, is related to the velocity of 
the light in free space c via 


a, ore (10) 
c 
The substitution of dB/dd from (9) into (6) yields the integral 
SO [1-vieBy-F (1-8) |Ear=o. (11) 
r| fa 2 p 


This integral was solved in Refs. (2) and (3) by assuming p constant and 
F = 2A(r/a)*, a power law, with a an arbitrary constant. To lift these 
restrictions and still solve (11) exactly, the following self-evident ex- 
pression is introduced: 
ry Op 
aan = p(roa, A) — p(ri, A) = 0. (12) 


ry 


This integral becomes less obvious and very useful when the derivative 
dp/or is performed with the help of (5), yielding 


re roF]| dr 
B-F-—|r—=0. 13 
J. | el p ( ) 


Combining (13) with (11), we arrive at a general expression 


ae aye f° 0-8) Sa 
T s 9 
ee eee (14) 


B r 
f ‘(145 r o) ear 
r1 2F or p 
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that like (11), relates the group delay t of a mode characterized by its 
mode parameter B (or propagation constant @), and its azimuthal mode 
number pv (hidden in p) to the profile function F and profile dispersion 
p. This expression is valid for any circular-symmetric fiber with isotropic 
dielectric and, in general, still cannot be solved exactly. However, if a 
particular family of fibers is considered that satisfies the condition 


pe 
2F or 


ue 
2 


= Dr), (15) 


D being an arbitrary function of \, the seemingly formidable right-hand 
side of (14) is reduced to 1/D and the group delay of the mode charac- 
terized by the mode parameter B is 


B 
1-— 
f= To. (16) 


These last two equations are the basic results of the paper. Equation 
(16) says that t, the group delay of a mode (or ray) is only a function of 
the mode parameter B and the dispersion parameter D. More important, 
the group delay is independent of the mode number (which means that 
modes with the same propagation constant have the same delay), it is 
independent of the profile function F, and it is independent of the profile 
dispersion p. Equation (16) is used in the following sections to study the 
impulse response of the fiber. 

On the other hand, eq. (15) is the recipe for the design of the fiber 
whose time response is given by (16). It can be solved in several ways 
depending on the control that the fiber designer has over F and its de- 
rivatives with respect to \. The least demands on these functions occur 
if the fiber is designed to operate at one wavelength only. Then, in (15), 
D becomes a constant, p is only a function of r, and the partial derivative 
of F is reduced to an ordinary one. Without introducing new symbols 
for D, F, and p, the simplified design formula is 


Tak 


2F dr 
cca 8; (17) 


8 
2 


This equation in turn can be solved in two ways. One way consists in 
prescribing the profile function F to satisfy, perhaps, requirements 
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different from modal dispersion. Then, the profile dispersion p must 
be tailored to satisfy (17). 

For example, assume the dielectric profile depicted in Fig. 1a. The 
profile function is . 


24 (—) for 0<r<a 
1 


F= (18) 


Q 


2A (-) ~ for ag<r<a, 
a 
where 
cal (cx —a ) 
Sear ) 12 (19) 
a . 
and the inequalities 
ay<a 
(20) 
ay >ag>O0 


guarantee that the profile looks indeed like that in Fig. la. 

Substituting (18) in (17) and assuming for the dispersion parameter 
D avalue Do that optimizes in some sense the impulse response of the 
fiber, the required profile dispersion turns out to be 


_2tay 


2 Do for r<dao 
p= a (21) 
ne for ag<r<a 
Do 


and is shown in Fig. 1b. 

This is an interesting example not only because it clearly illustrates 
the power of the theory even to design optimized fibers in which the 
profile function and dispersion are discontinuous, but also because it 
may be of practical interest. For example, by using an index-increasing 
dopant for r < ao and an index-decreasing dopant for r > ao the NA of 
the fiber can be increased, keeping at the same time its optimum modal 
dispersion. 

In the other way of solving (17), the profile dispersion p as a function 
F is assumed to be known, from experiment, and the index profile must 
be found from the integration of (17) that yields 


r=aexp fo a Os (22) 
F [2-—D(2—p)|F 


This result will be used in a more general way later, but if for the time 
being we prescribe p to be a constant Po, the profile function results: 
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(a) 


2—(2 + a7)/Do 


2—(2 + a4)/Do 


ARBITRARY 


ag a 
(b) 


Fig. 1—(a) Dielectric profile (solid line). (b) Profile dispersion. 


F=2A i) 
a 


a = D(2— Po) - 2. 


(23) 


(24) 


This last equation establishes the relation between the dispersion pa- 
rameter D of the fiber introduced in this paper and the a value so widely 
used in the literature?-? for fibers with constant profile dispersion Po. 
It follows from (17) that only if p is a constant, is the profile function F 
a power law (23). 
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The two solutions described require only the control over the profile 
function F and its first derivative with respect to A. But suppose that 
the fiber designer has also control over the second derivative. Then, to’ 
increase the range of wavelengths over which the fiber operates with low 
modal dispersion, he could demand, for example, that at a wavelength 
not only 


D= Do 
but also, as proposed by Kaminow and Presby,'® 


dD 
—— = 0. 25 
aX (25) 


This requires the simultaneous satisfaction of (17) and 
0 /r OF 
on (Far) 
dp 
an 


derived from (15) and (25). It is this last equation that implies the control 
over the second derivative of F with respect to X. 

It can easily be extrapolated that control over higher derivatives 
permits even further demands on D. In fact, if all the higher derivatives 
were controllable, D(A) could be chosen arbitrarily and the profile F 
would be the solution of the linear partial differential equation of first 
order (15) subject to the conditions of being zero at r = 0 and 2A(A) at 
r=a; the result iswell known!‘ from a mathematical point of view, but 
of limited importance from a practical point of view. 


= D9 (26) 


Ill. MINIMIZATION OF THE IMPULSE-RESPONSE WIDTH AND ITS 
SENSITIVITY TO ERRORS IN THE PROFILE 


The impulse- -response width is determined from (16) by finding the 
difference between flight times of the slowest and the fastest modes (or 
rays) for any given value of the dispersion parameter D. It is simple to 
find that the minimum time spread, Tmin, between those modes occurs 
if D is chosen 


Dp =14+V1— 2A. (27) 


In fact, the modes characterized by B = 0 and B = 2A are the slowest and 
arrive at the end of the fiber at 


bina = Dy (28) 
while the modes characterized by B =1—-V1— aa are the fastest and 


arrive at 
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2(1 — 2A)1/4 
tmin = . 
1+vV1—2A 


The time spread between them is the minimum impulse-response 
width 


(29) 


[1 -(1- 2QA)1/4]2 
Tmin = Emax — Emin = DS ee ee 
1+V1—2A 


Therefore, fibers with the same A have the same minimum impulse- 
response width tT, in, independently of their index profiles and profile 
dispersions, provided that they satisfy the design equation (17) with D 
substituted by the optimum value Do (27). 


(30) 


If 
A <1, 
(27) and (30) become 
Do x~2-—-A (31) 
A2 
Tmin = T = 0.61 A? us/km (32) 


for Ny = 1.46. 

To find the sensitivity of the impulse-response width to departures 
of the index profile from optimum, we calculate the ratio t/t min between 
the response width 7 for 


D = (1+6)Do, (33) 
where 
6<«1 


and the minimum response width 7,,in occurring for D = Do. After some 
straightforward calculations, 


2 
——=/l14 El (1+ V1—2A)? [1+ (1 - 2a)/a}2} (34) 
Tmin 
and for 
2 
= (1 + 1) (35) 
Tmin A 


It is known that the impulse-response width is indeed very sensitive to 
the choice of profile and more so for smaller A. If 6, the fractional de- 
parture of D from its optimum value, is equal to A, then the pulse width 
is nine times larger than 7,pin. 

The main results in the last two sections have been extended by Ar- 
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naud to optimize modal dispersion in fibers with noncircularly symmetric 
profiles. !5 


IV. THE RMS OF THE IMPULSE RESPONSE 


From the point of view of the maximum information-carrying capacity 
of a fiber, more significant than the impulse response width jis its rms 
width o,!6 since 1/4c is the repetition rate at which pulses can be trans- 
mitted with a reasonable loss penalty at the receiving end.! 

Let us calculate first the impulse response W(t) and then its rms width 
o assuming that: 


(i) The energy of the infinitely narrow impulse fed at one end of the 
fiber is equipartitioned among all modes. 
(11) All modes attenuate equally. 
(iii) The number of modes is so large that the discrete pulses arriving 
at the receiving end can be replaced by a continuum. 


The impulse response, then, is the rate of change of the number of 
modes reaching the end of the fiber, 


d Ymax 
Wit) = = f, uly, t)dy, (36) 


and its rms width is, by definition, 


i f f, ” W(t) W(te)(t1 — t2)2dt dt) 1” 


D) Ss 
f f, W(t 1)W(ts)dtidts 


To calculate W(t), the value of » given in (1) is substituted in (86) and 
the integration along » is carried through yielding 


o= 


(37) 


k2n? d "RB 
Wit ars B - F)rdr. 38 
(t) & Wes ( ) (38) 
The integral measures the energy arriving at the end of the fiber as a 
function of time and, since each contribution must be positive, the largest 
value that F' can reach is B. Therefore, the upper limit rg is the value of 
r that makes 


F(r, \) = B(A). (39) 


The explicit value of rg depends on the choice of fiber design. If the 
profile function F is prescribed, then rg is obtained by solving (39). If, 
on the other hand, the profile dispersion is prescribed, then 


2A dF : 
ra=aexp f° (2—-DQ-pIF (40) 
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follows from (22). 

Now the derivative with respect to time in (38) is carried out. The 
derivative of the integral is equal to the integral of the derivative since 
the terms that should contain the integrand times the derivatives of the 
limits are zero. Consequently, 


(knirp)? dB 
8 dt 


The reader interested in the explicit impulse response must substitute 
B in this expression with its time-dependent value obtained from 
(16). 

Replacing W(t) in (37) and also substituting the explicit value of t 
from (16), the rms width of the impulse response results in 


ff 1-x/D_1-y/D > ded 1/2 
T 0 7S Ve) ee 
Cee ee ee er 


V2 2A 
f f, rersdxdy 
0 


where x and y are dummy variables and r, and r, are given by (40) once 
B is substituted either by x or by y. It is easy to recognize in (42) that if 
A <«1andD ~ 2, the parenthesis is of the order of A? and c is propor- 
tional to A?T. 

Unlike the simple impulse-response width, the rms width o and the 
optimum value of D that minimizes it are dependent on the profile dis- 
persion p and the profile function F’. In general, the exact value of o and 
its minimizations must be found numerically, but we push the analysis 
a little further in the next section where some simplifying assumptions 
are made. 


W(t) = (41) 


, (42) 


V. APPROXIMATE RESULTS FOR RMS WIDTH OF THE IMPULSE 
RESPONSE, ITS MINIMIZATION, AND ITS SENSITIVITY TO PROFILE 
ERRORS 

Within the family of fibers described in the previous sections there 
is a large group of particular importance that encompasses many of the 
available fibers today. This group has small NA and its profile dispersion 
is almost constant with respect to r. To introduce these properties, we 
assume 


A<1, (43) 


then the profile dispersion is expanded in power series of the profile 
function F, 
foo F s 
=yP Ga) , 44 
p= LPs (on (44) 
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and since F is a function of r, the near invariance of p with r implies 


(55) «1. (45) 


u Ps\on 


s=l1 


5.1 Profile function 


Carrying this simplifying assumption to (22) by keeping only first 
powers of P,(s > 0), the profile function results: 


rasa)" free SL-G |p 4s 
where 
a = D(2— Po) — 2 (47) 


and D is still an arbitrary number. 
If the profile dispersion is constant, the summation in (46) disappears; 
then, and only then, will the profile function follow a pure power law. 


5.2 Minimization of the rms width of the impulse response and its sensitivity 
to profile errors 


We want to find omin, the minimum rms width of the impulse response 
possible, and D,, the optimum dispersion parameter for which omin is 
achieved. The optimum profile is obtained by substituting D with D, 
in (46). We are interested also in finding the sensitivity of o to small er- 
rors in the profile. 

To achieve these purposes o” is expanded in a power series about Dj, 
and only the first three terms.are kept, 


(D — Dy)? d20? 


do? 
2 = g2 i, + (D — Di) —~ + 
Tae 0 aD 2 dD? 


(48) 


The derivatives are to be taken at D = D,. Since by definition o? passes 
through a minimum of D = D,, the equation 


d 2, 
ai at D=D, (49) 


serves to determine the optimum dispersion parameter Dj. 
From (42), (48), and (49), we obtain with the help of (43) and (45) 
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A1l+ 2H 
D=2/ 1-5 Aat2)| (50) 
/ 
omin = T(AH)? Miata (51) 


(1+ 3H)(1 + 4H)(1 + 5H)1? 


o D-D,1+4H 2 
——= Se 1+ 3H)(1+ 5A), 52 
Omin v Gar eT, ( ) ( ) 
where 
-25P, (s?+s+6)H?2+ 8H+2 (53) 
=a “{I(s + 2)H + 1]? - Ak [(s + 3)H + 1]? — H7} 
and 
H=1-Pp». (54) 


The optimum value of the dispersion parameter D, is close to 2. The 
profile function that maximizes the information-carrying capacity of 
the fiber is obtained by substituting D with D, in (46). The dispersion- 
profile terms of order higher than zero appear in (50) only in the sum- 
mation > and are multiplied by A. Therefore, their contribution is small 
indeed and is neglected in (51) and (52). It is kept in (50) because, as will 
be seen later, small errors in the profile affect substantially the value of 
o. For > = 0, the substitution of (50) in (47) yields the same optimum 
a of Ref. 3. 

Consider now omin, the minimum rms width of the impulse response. 
From (51), we might be tempted to conclude that H = 1 — Po should be 
made small to decrease omin. However, the number of modes of the guide 
derived from (38) with the help of (46) is 








- (ey AH 


a 50 
2 1+ ea! 


Therefore, if the number of modes of the fiber is to be kept constant, omin 
can be decreased by making 


(1 + H)5/2 
(1 + 3H)(1 + 4H)(1 + 5H)! 


small and this is achieved by choosing H large, not small. 

The following table contains the minimum rms per unit length of fiber, 
Omin/L, and the concomitant pulse repetition rate L/4omin for different 
values of H as derived from (58), assuming N, = 1.46. 
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L 





H GainlL PRR = 
Omin 
us/km Mb- km/s 
1 0.14A2 1.79/A2 
2 0.16A2 1.56/A2 
oo 0.18A2 1.38/A2 


For H > 1, the pulse repetition rate is fairly insensitive to the value of 
H.¥or H = 1 and A = 0.01, the pulse repetition rate is ~18 Gb/km/s. This 
information-carrying capacity is only 33 percent smaller than that of the 
“ideal profile” reported by Cook.!” 

Let us turn now to the sensitivity of the rms width to errors in profile 
(52). Again, for H > 1, this result is insensitive to the value of H; in- 
deed 
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(56) 





For H = Land (D — D,)/Dy,, the fractional departure of D with respect 
to the optimum D, equal to A, the rms width o is about 8.4 times wider 
than omin. As in the case of the pulse width, the rms width is very sensi- 
tive to profile errors. 

A fiber designed to minimize the rms width (D = Dj) has only 30 
percent more information-carrying capacity than a fiber with the same 
A designed to minimize the impulse-response width (D = Do). 


VI. CONCLUSIONS 


For a vast class of circular-symmetric fibers made of isotropic di- 
electrics, simple and fundamental design criteria that minimize the 
impulse-response width due to modal dispersion at one wavelength have 
been found. This minimum width (30) is only a function of the NA and 
the time of flight along the axis. Therefore, if properly designed, a fiber 
with arbitrary profile dispersion has the same minimum impulse-re- 
sponse width as another fiber with the same NA and no profile dispersion. 
Their information-carrying capacity is about 1.4/A? Mb - km/s. The fiber 
engineer has a substantial freedom of choice to reach that optimum 
design: the profile dispersion may be arbitrarily chosen but then the 
index profile is uniquely determined by (22); or symmetrically, the index 
profile may be arbitrarily chosen and then the profile dispersion must 
satisfy (17). Only if the profile dispersion is a constant does the optimum 
dielectric profile that minimizes the impulse response follow a power 
law. 
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The profile dispersion entails the first derivative of the index with 
respect to the wavelength. If the second derivative can be controlled, then 
the minimization of the impulse-response width can be achieved at two 
neighboring wavelengths. This broadbanding of the fiber response can 
be expanded even further if higher derivatives are under control. 

The width of the impulse response is very sensitive to errors in the 
fiber design. A fractional error of A in the dispersion parameter of the 
fiber makes the response about nine times wider than the minimum as 
seen from (35). 

Only a marginal increase of about 25 percent in the information- 
carrying capacity of the fiber is achieved if the rms width of the impulse 
response is minimized instead of minimizing the impulse-response 


width. 
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An Experimental Optical-Fiber Link for Low-Bit- 
Rate Applications 


By W. M. MUSKA 
(Manuscript received July 15, 1976) 


The design, construction, and performance of a practical optical link 
with interface circuits for coding and decoding 1.5-Mb/s, bipolar, digital 
signals are described. The optical devices used are light-emitting diodes 
and PIN photodiodes. A feedback or “transimpedance”’ preamplifier 
that incorporates a silicon junction-field-effect transistor is used in 
the receiver, which has a sensitivity of —57.2 dBm average optical power 
for a 10~° bit-error-rate. The receiver demonstrates an optical power 
dynamic range of about 28 dB without requiring automatic gain 
control. Timing recovery is accomplished by a simple, conventional 
technique. 


I. INTRODUCTION 


Experimental repeaters for regenerating digital signals from a few 
Mb/s to a few hundred Mb/s for optical-fiber transmission have been 
reported.'~’ In this paper, we report the results of an experiment in which 
a practical optical repeater of simple design and high performance with 
interface circuits for coding and decoding 1.5-Mb/s bipolar signals was 
constructed and evaluated. Such a low bit rate, bipolar, digital format 
-(DS1) is presently used in telephone systems for the transmission of 
multiplexed, digitally encoded voice signals over copper twisted pairs. 
The optical system was designed to be transparent to the bipolar format 
of these electrical signals. It utilizes a directly modulated light-emit- 
ting-diode (LED) source and a PIN photodiode, incorporates a simple 
timing recovery circuit, and demonstrates receiver sensitivity that ap- 
proaches theory. 

The bipolar format consists of “zeros” and 50-percent duty cycle al- 
ternating positive and negatives “ones”; it, therefore, has three levels 
with a zero dc component. A straightforward scheme to translate the 
bipolar format to a unipolar format without increasing the bit rate might 
involve simple full-wave rectification. However, this is unacceptable in 
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practice because violation of the bipolar format is used in system 
maintenance functions as an indication of system performance. 

If the bipolar signal is transmitted as a three-level optical signal rather 
than a binary signal at the bipolar bit rate, 3 dB more average optical 
power would be required at the receiver detector for a given error-rate 
(assuming a thermal-noise dominated receiver). Implementing the 
three-level system would be more difficult, due to the nonlinearity of 
the LED. The three-level system also requires separate regenerators in 
the receiver for +1 and —1 pulse regeneration. If the bipolar signal is 
coded in a two-bit-for-one-bit manner (twice the bit rate), the receiver 
sensitivity would suffer a degradation of about 4.5 dB (PIN photodiode 
with FET amplifier input device)? compared to a binary format system 
at the bipolar bit rate. Although this system is about 1.5 dB less sensitive 
than the three-level system, repeater circuitry is simplified since the 
effects of LED nonlinearity are not important and only one pulse re- 
generator is required. 

More efficient coding schemes are possible (e.g., 3-bit for 2-bit)? that 
allow the bipolar information to be transmitted at less than twice the 
information bit rate; however, the circuitry required to implement these 
schemes is substantially more complicated. 

The simple coding format?!” used in our experiment is shown in Fig. 
1. A positive one is coded into two consecutive 3-Mb/s ones (non-re- 
turn-to-zero in this case) within a coding frame, a zero is coded as a one 
followed by a zero within a coding frame (this assignment is arbitrary; 
it could be 0, 1, instead), and a —1 is coded as two consecutive zeros. This 
coding format maintains desirable features of the bipolar signal, such 
as dc balance and substantial redundancy for error correction or moni- 
toring. In addition, the code allows a receiver that is ac-coupled to have 
a constant threshold level at zero volts for a large range in received power 
level without automatic gain control. The pulse-transition density of the 
coded signal is very high, thus allowing a timing-recovery circuit in the 
receiver that is less critical than those encountered in more conventional 
binary systems. 

The optical system is designed around an LED light source and a PIN 
photodiode. These devices are expected to be less expensive and require 
less control circuitry than lasers or avalanche detectors. 


Il. DESCRIPTION 


A unidirectional optical-fiber link is shown in Fig. 2. The coder- 
transmitter converts the 1.5-Mb/s incoming bipolar signal into a 3-Mb/s 
unipolar optical signal for transmission via optical fibers. At the receiving 
terminal, the optical signal is regenerated and decoded back to the 
1.5-Mb/s bipolar format. Depending upon the distance between termi- 
nals, one or more line repeaters may be needed. A line repeater might 
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Fig. 1—Coding format. 


include the receiver and driver shown in Fig. 2. and a simple regenerator 
since there is no need to decode or encode the unipolar signal. 


2.1 Coder-transmitter 


The circuits between the incoming electrical bipolar signal and the 
optical fiber consist of a bipolar-to-unipolar converter, a driver for the 
LED, and the LED. The bipolar signal is coupled to the converter by a 
transformer with a center-tapped secondary. Each half of the secondary 
is half-wave rectified; one rectifier generates a pulse when a “plus-one” 
appears at the input and the other generates a pulse when a “minus-one”’ 
appears at the input. 

These outputs are then encoded into the 3-Mb/s format by TTL cir- 
cuitry realized in two dual-in-line packages. The encoding process re- 
quires a clock signal, which is obtained from a timing-recovery circuit 
similar to that used in conventional digital repeaters.!! 

The digital output voltage of the converter is converted into 220-mA 
peak current pulses, which directly modulate the Burrus-type, diffused 
junction LED. A non-return-to-zero format is used since these LEDs are 
basically peak-power limited and the tolerable repeater span increases 
with average power when fiber-delay distortion is negligible. The 
bandwidth of these devices is more than adequate for this bit rate. 
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Fig. 2—Electronics of the optical-fiber link. 


2.2 Receiver-decoder 


The receiver and decoder-regenerator converts the unipolar signal 
on the optical fiber to a bipolar electrical signal. The input fiber is cou- 
pled to a PIN photodiode, which converts the optical pulses to electrical 
pulses which are then amplified sufficiently to be processed. The sen- 
sitivity of the receiver is determined by the thermal and shot noise 
generated in the “front-end” amplifier and the quantum efficiency and 
shot noise of the PIN photodiode. It has been shown that the total noise 
of the amplifier and diode may be referred to the input of the amplifier 
as a shunt current generator.!? In the receiver described, the dominant 
source of noise is the silicon junction-field-effect transistor (Si JFET) 
input device since the leakage-current shot noise for most PIN photo- 
diodes with active-area diameters less than 1.25 mm is small enough to 
be neglected. The mean-square current of the effective input noise 
current generator has been shown to be! 
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where 


k = Boltzmann’s constant 
T = absolute temperature 
fy = bit rate = 3 MHz 
£m = transconductance of the FET ~ 5 mA/V 
Cr = Caiode + Cstray + Cin FET = 8 pF 
R, = parallel equivalent resistance at the input to the FET 
(diode bias return, FET bias, etc.) 
Yl ~ 0.7 for typical Si JFETs. 


The parameters I and J, are weighting factors, which are determined 
by optical input and amplifier output pulse shapes;!? in this case the 
input pulse is NRZ and rectangular and the output pulse NRZ and raised 
cosine; therefore, 


Io = 0.55 
Io = 0.085. 


If the amplifier input device is an FET, it is possible to choose a value 
for R,, such that the term associated with Io is negligible compared to 
the term associated with I», in which case an FET should be chosen to 
maximize Vg,,/Cr. 

The above analysis is for a non-feedback amplifier, which has an input 
time constant much larger than the signal period, 1/fy. In the present 
work, a “transimpedance” or feedback amplifier is used. Here, the output 
voltage is 





where 


I, = photodiode signal current 

Zt = feedback impedance = 1/(1/R; + jwC;) 

Cy = total shunting capacitance across R; 

R; = feedback resistance 

A = open-loop amplifier voltage gain 

B= Zin/(Zin + Z,) 
Zin = total impedance at the input of the open-loop 

amplifier 
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It can be shown that the feedback resistor R; plays the same role as R, 
in relation to noise. 

Since, in this case, C; is about (1/80)C;, the time constant R;C; is 1/80 
that of R,C; in the non-feedback design, and the signal is integrated over 
far fewer signal periods. There must be sufficient linear gain before pulse 
shape restoration, so that the amplifier following equalization does not 
affect the signal-to-noise ratio. In the case of the non-feedback amplifier, 
the integrated peak output voltage excursion may be many times that 
for a single pulse, in which case the optical power dynamic range of the 
amplifier is limited. The transimpedance amplifier, on the other hand, 
does not experience this problem, thus allowing a substantially greater 
dynamic range. 

Following the pulse equalizer, the signal is further amplified by an 
amplifier that is linear at very low input power levels but acts as a sym- 
metrical limiter at higher received optical power levels without degrading 
the performance. This feature eliminates the need for an automatic gain 
control. 

The signal is bandlimited by a third-order Butterworth filter to reduce 
high-frequency, out-of-band noise, and then threshold detected by a 
zero-crossing detector. 

Due to the coding format, the dc component of the received electrical 
signal voltage is always one-half the peak, thus the proper threshold level 
at the threshold detector is always zero when ac coupling is used at the 
input of the threshold detector. 

The signal at the output of the zero-crossing detector, which has been 
quantized into two discrete voltage levels, has timing-jitter due to the 
presence of noise. This jitter is removed, as shown in Fig. 3, by sampling 
the signal at times which are determined by a 1.5-MHz, low-noise clock 
signal, which is recovered from the received signal by a high-Q, paral- 
lel-resonant circuit. This circuit is sustained in oscillation by 3.3-mi- 
crosecond pulses, which are generated by a monostable multivibrator 
whenever a zero-one transition occurs in the received signal (Fig. 3, 
waveforms A-C). Since a zero-one sequence in the encoded signal only 
occurs at the beginning of a coding frame, the clock is synchronized with 
the frames. The 3-Mb/s received signal is demultiplexed into two 1.5- 
Mb/s signals by two sample and store circuits, one of which samples the 
first bit in a coding frame on positive-going clock transistions (Fig. 3, 
waveform D), and the other samples the second bit in a frame on nega- 
tive-going clock transistions (Fig. 3, waveform E). These two 1.5-Mb/s 
signals are then reconstituted into the original DS1 (1.5-Mb/s, 50-percent 
duty-cycle, bipolar) signal by comparing the states of the two demulti- 
plexed signals at the end of each coding frame (Fig. 3, waveform F). Pulse 
regeneration and decoding to bipolar is accomplished with TTL circuits 
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Fig. 3—Decoding process. 


consisting of 3 D-type flip-flops, 3-2 input NAND gates, and an inverter. 
The output of the decoder is coupled to a twisted pair line. 


ll. SYSTEM PERFORMANCE 
3.1 Evaluation procedure 

To evaluate the system error performance, dynamic range, and timing 
jitter, the coder-transmitter and the receiver-decoder were physically 
separated and optically coupled through an air path into which the re- 


quired amount of attenuation was placed. The two sections were then 
packaged in a plug-in module of dimensions, 1% inches X 8'4 inches X 
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9 inches, and optically coupled via an optical-fiber path into which a 
controlled amount of attenuation was placed. No degradation in per- 
formance was observed. 


3.2 Error performance 


The theoretical average optical power required for a given error 
probability is! 


hy 
P ac PTEN Qn, 1/2, 
ne 


where 


P = average optical power required to produce an error 
probability of P, 

h = Planck’s constant 

v = optical frequency 

n = photodiode quantum efficiency 

e = electron charge 

Q = V2 erfc—! (2P,) 

P. = probability of error. 


The photodiode used in the error probability measurements had a 
quantum efficiency of about 65 percent. 

The theoretical and experimental error probability curves are pre- 
sented in Fig. 4. The 0.5-dB discrepancy is due to noise contributions 
of amplifier stages following the FET input stage and measurement er- 
rors. The data point at the 1.4 X 107° rate was obtained by counting the 
errors that occurred in 12 time periods of 10° time slots each (~5 min.). 
During one of these periods, a burst of 30 errors occurred due to elec- 
tromagnetic interference and the data obtained during this period were 
omitted. 

The curves in Fig. 4 indicate the error probability at the coded 3-Mb/s 
rate, the error probability at the 1.5-Mb/s rate would be about a factor 
of two higher (some error correction takes place in the decoder); however, 
this increase could be compensated for with about a 0.2 dB increase in 
optical power. 


3.3 Dynamic range 
Optical power dynamic range is described as 
Pa (dB) a Prax (dB) _ Prin (dB), 


where 
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Fig. 4—-Experimental and theoretical error vs. received average optical power. 


Pmax = Maximum received optical power level before nonlinear 
effects degrade performance. 

Pin = optical power required to maintain a 107° error 
probability. 


The optical power dynamic range for the receiver described here is about 
28 dB. 
3.4 Timing jitter 

The phase jitter in the decoded pulse stream was observed to be about 
10 degrees RMS at a received optical power level of about —57.8 dBm. 
3.5 Waveforms 


The oscilloscope pictures in Fig. 5(a—e) show waveforms at several 
points in the system at various input power levels. The ‘“‘eye diagram” 
in Fig. 5(a) is taken at the output of the “front-end” amplifier and il- 
lustrates the signal integration due to the capacitance across the feedback 
resistor. Fig. 5(b) and (c) were taken at the input to the filter, Fig. 5(b) 
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Fig. 5—Receiver waveforms (voltages: 0.1 < actual value) 


was taken at a received optical power level of —47 dBm, and Fig. 5(c) was 
taken at an optical power level of —28 dBm when the second amplifier 
acts as a limiter. 

Figure 5(d) demonstrates the “eye,” after filtering, at a —57.2 dBm 
optical power level. The waveforms from top to bottom in Fig. 5(e) are: 
the bipolar signal at the input to the coder-transmitter, the 3-Mb/s re- 
ceived signal at the output of the filter, and the regenerated and decoded 
bipolar signal. 


Power consumption 


Total power consumption for the transmitter and receiver was about 
1.75 watts. 
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IV. SUMMARY 


A practical optical-fiber link of simple design for faithful transmission 
of 1.5-Mb/s bipolar, digital signals has been built and tested. High per- 
formance is achieved through the use of a transimpedance preamplifier 
that affords receiver sensitivity that approaches the theoretical limits 
imposed by commercially available silicon junction-field-effect tran- 
sistors. 
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Propagation of High-Frequency Elastic Surface 
Waves Along Cylinders With Various Cross- 
Sectional Shapes 


By J. A. MORRISON, J. B. SEERY, and L. O. WILSON 
(Manuscript received March 24, 1976) 


Elastic surface waves, or Rayleigh waves, are disturbances that travel 
over the stress-free surface of an elastic solid, and whose amplitudes 
decay rapidly with depth into the solid. Earlier mathematical results. 
are used to study numerically the properties of these waves on specific 
cylindrical objects that might be used as acoustic topographic wave- 
guides. The lowest-order mode is investigated for cylinders with strictly 
nonconstant curvature. Mode confinement and its dependence on such 
things as cylinder shape and the value of the frequency parameter are 
studied. Phase and group velocities are also computed. Mode behavior 
1s studied in the transition region between the case of cross-sectional 
boundary curves of nonconstant (and not “almost” constant) curvature, 
for which the modes are localized, and the case of constant curvature, 
for which they are not localized. Some higher-order modes are inves- 
tigated for the rounded wedge. 


I. INTRODUCTION 


Elastic surface waves, or Rayleigh waves, are disturbances that travel 
over the stress-free surface of an elastic solid, and whose amplitudes 
decay rapidly with depth into the solid. In a series of earlier papers,!~° 
we developed and applied some mathematical techniques to describe 
the propagation of high-frequency elastic surface waves along cylinders 
of general cross section. Our intent was to learn more about the prop- 
erties of such waves traveling down cylindrical objects that might be used 
as acoustic topographic waveguides. In this paper, we use our earlier 
mathematical results to study numerically the properties of elastic 
surface waves on certain specific cylindrical objects of interest. We treat 
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cylinders roughly corresponding to an elliptical bore, an elliptical rod, 
a wedge with a rounded tip, and a flat plane with a rounded ridge on it. 
The elastic medium is assumed to be homogeneous and isotropic. 

The earlier papers discussed two approximate high-frequency de- 
scriptions of the surface-wave behavior: an asymptotic approximation 
and one which we termed a surface-wave approximation. The analysis 
involved a scalar wave equation, a vector wave equation, and rather 
complicated boundary conditions. Since the analysis was cumbersome, 
a simpler scalar “model problem” was first investigated by Morrison.! 
The techniques he developed had counterparts in the full elastic prob- 
lem, which was treated by Wilson and Morrison? in the high-frequency 
asymptotic approximation, designated by A (as depicted in Fig. 1). The 
lowest-order surface-wave mode was investigated in almost as much 
detail as that for the scalar problem, but because of the algebraic com- 
plexities, the higher-order modes were less completely analyzed. 

For the scalar problem,! Morrison had also obtained a surface-wave 
approximation describing the high-frequency behavior of the surface- 
wave modes. He then derived the analogous approximate equations for 
the high-frequency behavior of the elastic surface-wave modes.® Unlike 
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Fig. 1—High-frequency approximation chart. 
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the asymptotic approximation, A in Fig. 1, the surface-wave approxi- 
mation, designated by B in Fig. 1, was capable of describing the mode 
behavior at high frequencies in the transition region between the case 
of cross-sectional boundary curves of nonconstant (and not “almost”’ 
constant) curvature, for which the modes are localized, and the case of 
constant curvature, for which they are not localized. (The asymptotic 
approximations had required boundary curves with strictly nonconstant 
curvature.) Also, this description gave a more complete analysis of the 
higher-order modes. 

In Section II, we exhibit the cross-sectional boundary curvature 
functions used in our numerical investigations and explain why we chose 
those specific functions. 

Section III is devoted to a numerical treatment of the high-frequency 
asymptotic results. The lowest-order mode is investigated for cylinders 
with strictly nonconstant curvature. We learn about the phenomenon 
of mode confinement, and its dependence upon such things as the shape 
of the cylinder and the value of the high-frequency parameter x used 
in the asymptotic expansions, A in Fig. 1. We also compute the phase 
and group velocities. 

In Section IV we make a similar investigation using the surface-wave 
approximation, B in Fig. 1. When possible, the results are compared with 
the asymptotic results. Particular attention is paid to the mode behavior 
in the transition region described earlier, and to the behavior of 
higher-order modes. We also compare our results with exact theoretical 
results for the circular bore.4 

In Section V, we summarize our findings. 


ll. THE BOUNDARY CURVES 


We wished to investigate numerically the properties of disturbances 
propagating along the surfaces of various cylindrical objects. The mo- 
tivation for our particular choices of cross-sectional boundary curves 
came from our earlier high-frequency asymptotic results,” which could 
be applied to a cylinder with an open boundary curve whose curvature 
attains its algebraic maximum at a single point, and which could also be 
applied to a cylinder that has a closed boundary curve which is symmetric 
and whose curvature attains its algebraic maximum at two points. We 
decided to consider disturbances propagating along objects roughly 
corresponding to an elliptical bore, an elliptical rod, a wedge with a 
rounded tip, and a flat plane with a rounded ridge on it. 

The exact forms of the chosen boundary curvature functions were 
suggested by the analytical form of the displacement function obtained 
in our high-frequency asymptotic results.* The high-frequency behavior 
of the disturbance can be determined in the vicinity of the cylinder 
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surface. As is shown in Appendix A, the disturbance corresponding to 
the zeroth-order mode, when evaluated at the surface of the cylinder, 
can be expressed as 


c e'(bz—wt) A 
“pb OD) BEM Z-0 
b2 2 C 
oy re Clea Wea . C(n) 
[3 2b? IQ Tene + ib | 1+ 5 a |). (1) 


Here u is the displacement, @ is the propagation constant, z measures 
distance along the generators of the cylinder, w is the frequency, and t 
is the time. A right-handed coordinate system is used, with unit vectors 
n, t, and k in the directions of the inward normal, tangent to the cross- 
sectional boundary curve, and along the generators of the cylinder, re- 
spectively. Here = = n/€, where n represents distance from the surface 
along the inward normal and is a characteristic length; also n = s/@, 
where s is signed arc length along the boundary curve. The normalized - 
unit of length, corresponding to the characteristic length @, is depicted 
in Figs. 2, 3, and 4 for the particular cross-sectional boundary curves 
considered. The quantities b, a7, az, and P are constants defined in 
Appendix A, and })(0) is a normalization constant. We have 


x =wl/er > 1, (2) 


where cr is a constant representing the transverse wave velocity of the 
medium. Thus, the parameter x is proportional to the frequency w and 
is assumed to be large. The functions C(n), F(), G(m), and J(7) all involve 
the curvature function K(n) = ¢x(s); as is shown in Appendix A, C(n), 
F(n), and G(n) are defined as integrals of certain functions of the cur- 
vature function. It was to evaluate these integrals and the integral of the 
curvature function analytically that we chose the specific curvature 
functions. 

The boundary curves are then given by the functions X (7) and Y(n) 
defined by® 


oS cos §. K(ods, 
dn 0 
iy 
oF an f, "K(odé. (3) 
dy 0 
To describe an ellipse-like bore, we set 
Ki(y)=—-1+2kcos2n, OSk 2%, (4) 
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Fig. 2—Boundary curves corresponding to the ellipse-like bore and rod for various values 
of k. 
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Fig. 3—Boundary curves corresponding to the rounded wedge for various values of 
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Fig. 4—Boundary curves corresponding to the ridged plane for various values of k. 


while 
Ko(n) = 1 + 2k cos 2n, Osksk (5) 


corresponds to an ellipse-like rod. Equations (4) and (5) actually describe 
the same boundary curve with circumference 27. Notice that in both 
instances 7 = 0 corresponds to a point of maximum algebraic curvature. 

Also, the curvature is negative for the bore, because of our convention 

that n is directed into the region. The curvature function 


K3(n) = k sech?n, Osks7/2 (6) 


corresponds to a wedge-like object similar to a hyperbolic cylinder. Fi- 
nally, 


K4a(n) = k[e? — tanh?n(sech?n + €)2], e= l, (1 + Vz) (7) 


describes a rounded ridge on a plane. The value for ¢ is found from the 
condition {9 K(n)dy = 0. As with the other cases, there is also a re- 
striction on k. It is complicated so we do not write it here. In each case, 
k is a parameter that can be varied. 

The boundary curves corresponding to the curvature functions in (4) 
to (7) were obtained numerically. In Fig. 2, we show the boundary curves 
corresponding to the ellipse-like bore and rod for various values of k. The 
unit of length is indicated in the figure. Notice that k = 0 corresponds 
to a circular bore or rod. The dot represents the point 7 = 0 for the bore 
corresponding to k = 0.1, indicating that the points of maximum cur- 
vature lie on the Yj axis. For the bore, the vectors n, and t, indicate the 
directions of the normal into the region and the tangent to the curve. 
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Similarly, the circle represents the point 7 = 0 for the rod corresponding 
to k = 0.5, indicating that the points of maximum curvature lie on the 
Y> axis. Note that the axes for the rod have been rotated by 90 degrees 
from those for the bore. For the rod, the vectors ng and tz indicate the 
directions of the inward normal and the tangent to the curve. In Figs. 
3 and 4, we show the boundary curves for the rounded wedge and ridged 
plane, respectively, for various values of k. The unit of length is indicated, 
as are the unit vectors n and t. 

At high frequencies, the asymptotic result (1) does not hold in the 
transition region between the case of cross-sectional boundary curves 
of nonconstant (and not “almost” constant) curvature, for which the 
modes are localized, and the case of constant curvature, for which they 
are not localized. However, a refined surface-wave approximation 
equation, to be discussed in Section IV, does give results in this transition 
region. 


lll, ASYMPTOTIC RESULTS 


In this section, we present high-frequency asymptotic results based 
on evaluation of (1) for the bore, rod, wedge, and plane with a ridge. In 
each case, the results are for the fundamental, or zeroth-order, mode. 
The bore and the rod have closed boundary curves which are symmetric, 
and for which the curvature attains its algebraic maximum at two points. 
For such cylinders, the expansion (1) corresponds to two modes, the 
zeroth-order symmetric one and the zeroth-order antisymmetric one, 
for which the values of the propagation constant @ differ by only an ex- 
ponentially small amount.!:? This expansion is about the point of max- 
imum algebraic curvature at n = 0 and is valid for |n| < 7/2. There is an 
analogous expansion about the point of maximum algebraic curvature 
at » = a which is valid for |y — z| < 2/2. Each expansion is not expected 
to be precise in regions where the disturbance is very small and the two 
modes differ. It is necessary to ensure that the disturbance is confined 
to regions near points of maximum algebraic curvature so that it is indeed 
small where the two modes are known to differ. From (1), (29), (33), and 
(35), this can be viewed as a requirement that the frequency parameter 
x be sufficiently large and that the deviation of the curvature from a 
constant value not be small. 

Equation (1), which is valid on the surface of the cylinder, was ob- 
tained from (39) in Appendix A, which holds also near the surface. In 
the derivation of the latter equation, it was necessary to assume that if 
the center of curvature for a point on the cross-sectional boundary curve 
lies within the region defining the cylinder, then the disturbance must 
be negligible at that point. This means that for the rod, wedge, and ridged 
plane, the results are applicable only if the frequency is high enough that 
the disturbance is confined close to the surface of the cylinder. Mathe- 
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matically, this means that each exponential term in (39) must be very 
small when evaluated at the value of % corresponding to the minimum 
radius of curvature. These conditions result in an approximate lower 
bound on the frequency parameter x, namely that 


x Z 10 K(0)/ar. (8) 


The constant a7 depends upon the Poisson ratio o. For values of o be- 
tween 0 and 44, it turns out that ar ranges between 0.56 and 0.31, re- 
spectively. Most of our numerical results were obtained with o = 0.16974, 
which was taken as the Poisson ratio for fused silica.® For this value of 
a, we have ar = 0.47. 

As can be seen from (1), the n and k components of the ‘displacement 
Un(0,n) and uz, (0,n) each contain a factor [1 + 44C(n)(Px)—!], whereas 
the t component u;(0,n) contains a factor I(n)(P/x)!/2. Then, in the 
lowest-order asymptotic approximation, A.0 in Fig. 1, the t component 
of the displacement does not even appear. To this order, the solution is 
like that for Rayleigh waves traveling on the surface of a plane infinite 
half space except that it is multiplied by a factor that describes the 
confinement of the disturbance due to the cylinder curvature. In the 
next-order asymptotic approximation, A.1 in Fig. 1, the effect is to 
multiply this solution by an additional factor and to add a t component 
of displacement. Since it turns out computationally that u; is a few 
percent of the size of u,, or uz, we shall concentrate our attention on u,, 
and uz. Since these two components are proportional to each other, it 
suffices to treat the quantities 


So(n) = F(n) exp[—(Px)1/2G ()], (9) 
Si(m) = So(n)[1 + %C(n)(Px)-¥7], (10) 


which are proportional to the n and k displacement components in the 
lowest-order and next-higher-order asymptotic approximations, A.0 and 
A.1, respectively, and which are normalized to unity at 7 = 0. We shall 
mostly discuss the more accurate approximation S,(n). 

Figures 5 and 6 illustrate some results obtained for the ellipse-like bore 
whose cross-sectional curvature function is given by (4). Because of the 
symmetries of the boundary curve, it is only necessary to consider values 
of n between 0, which corresponds to a point of maximum algebraic 
curvature, and 7/2, which corresponds to a point one quarter of the way 
around the curve. In both figures, we plot the first-order asymptotic 
approximation S;(n) from 7 = 0 to a value of 7 less than 2/2 for which 
the disturbance is relatively small. 

In Fig. 5, we show S;(n) for several values of the frequency parameter 
x. The constant k in (4) was chosen to be 0.5; the Poisson ratio was chosen 
to be 0.16974, corresponding to fused silica. If we take the transverse wave 
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Fig. 5—First-order asymptotic approximation S;(7) for the bore, with k = 0.5, plotted 
as a function of y for various values of the frequency parameter y and Poisson ratio o = 
0.16974. 


velocity for fused silica to be 3764 m/s, then for 2 = 3 X 10-4 m the fre- 
quency corresponding to x = 40 is approximately 80 MHz. It is strikingly 
apparent that the disturbance is indeed confined to a region near the 
point of maximum algebraic curvature y = 0. There is, of course, similar 
confinement to the region near y = z. Such confinement near a point of 
maximum algebraic curvature shows up in a similar manner in the 
computations for rods, wedges, and ridges on planes. We also see that 
the confinement becomes even more pronounced as the frequency pa- 
rameter y increases. 
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Fig. 6—First-order asymptotic approximation S(y) for the bore, with k = 0.5, plotted 
as a function of y for frequency parameter x = 80 and various values of the Poisson ratio 
oO. 
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As the ellipse-like bore becomes more like a circular bore (k decreases), 
the confinement near points of maximum algebraic curvature decreases. 
This is to be expected since in the limiting case of a perfectly circular 
bore, there is no such confinement at all. We do not include a figure here 
to exhibit the way the confinement changes with k, as the effect will be 
vividly demonstrated later in Section IV where surface waves on “al- 
most” circular bores are discussed. 

In Fig. 6, we fix y = 80 and k = 0.5, and show S,(n) for various values 
of the Poisson ratio c. The confinement increases as o decreases. 

These results are quite representative of all those we obtained for the 
bore, rod, wedge, and plane with a ridge. Other curves for S;(n) are 
qualitatively very similar. For example, in Fig. 7, we show S;(n) for a 
plane with a ridge on it. In this case, the curvature function is given by 
(7), with k = 3. We fixed o = 0.16974 and varied the frequency parameter 
x. Because other results are so similar, we do not show any specific curves 
for the rod or wedge. 

We next compare the lowest-order asymptotic approximation So(7) 
with the first-order asymptotic approximation S;(n). Here, some dis- 
tinctions do arise in our computations for the various cylinders. In 
considering the ellipse-like bore, we find that for values of the parameters 
in the ranges previously discussed, the curves So(7) are hardly distin- 
guishable from the curves S;(n). This is not always the case for the rod, 
wedge, and ridge on a plane. To give an example for which the first-order 
correction term is significant, we show in Fig. 8 the functions So(n) and 
S1(m) corresponding to a rounded wedge, with x = 40, k = 1.0, and o = 
0.16974. Similar results can be obtained for the rod and the plane with 
a ridge on it when the frequency parameter x is in the lower part of the 
range being considered. In all cases, the first-order correction becomes 
noticeably smaller as x is increased. We expect the second-order cor- 
rection to be negligible. 

The normalized phase and group velocities are 


__o _ (,,,a48\7} 
a ae (era.) ; at) 


where cr is the transverse-wave velocity. From the asymptotic results,” 
it is found that 
1 doP (—doP)/2 (dP? — 2b2Q?) 


a sire 12 
Wp b  2b3x 2633/2 Ab5y2 (12) 
and 

_1 (-deP)? | @Q 
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Fig. 7—First-order asymptotic approximation S;(y) for the ridged plane, with k = 3, 
plotted as a function of 7 for various values of the frequency parameter x and Poisson ratio 
o = 0.16974. 


where the constants b, do, do, P, and Q are defined in Appendix A. In 
particular, 1/b = cr/cr = wp is the normalized Rayleigh wave velocity, 
and wp is the solution of equation (24) which satisfies 0 < wr < 1. Both 
the phase and the group velocity asymptotically approach the Rayleigh 
wave velocity and, consequently, we define the normalized differential 
phase and group velocities by 


bp = Wp — 1/b, bg = Wg — 1/b. (14) 

The asymptotic approximations 6p1, 5p2, and 6p3 to dp are obtained 

by retaining terms through orders x~!, x~3/2, and y~? respectively in the 
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Fig. 8—Comparison of lowest- and first-order asymptotic approximations So(y) and 


S1(n) for the rounded wedge, with k = 1.0, plotted as functions of 7 for frequency parameter 
x = 40 and Poisson ratio o = 0.16974. 
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expansion (12). Similarly, the approximations 6,2 and d,3 to dg are ob- 
tained by retaining terms through orders y~?/2 and y~? respectively in 
the expansion (13). 

In Fig. 9, we show the approximations 6,2 and 673 to the normalized 
differential phase velocity as a function of the frequency parameter x 
for the ellipse-like bore, with k = 0.5 and o = 0.16974. We do not plot 6p, 
as it can be shown to be identically equal to zero for this case. Notice that 
the convergence is quite good. This is also true for the approximations 
6,2 and 6,3 to the normalized differential group velocity, which are shown 
in Fig. 10. 

For cylinders of other cross-sectional shapes, the convergence is not 
always so good, particularly for the differential group velocities. In Tables 
I(a) and I(b), we show 51, dpa, 5p3, and dg, dg3, respectively, for the bore, 
rod, wedge, and ridged plane; here x = 80, o = 0.16974, and k varies. The 
convergence improves as x increases. 

In Figs. 11 and 12, we show one additional set of approximations to 
the differential phase and group velocities, respectively, as a function 
of x. Here, the curves are for the wedge, with k = 1.0 and o = 0.16974. 
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Fig. 9—Asymptotic approximations 6p2 and 6p3 to the normalized differential phase 
velocity as a function of the frequency parameter x for the bore; k = 0.5 and Poisson ratio 
o = 0.16974. 
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Fig. 10—Asymptotic approximations 6,2 and 6,3 to the normalized differential group 
velocity as a function of the frequency parameter x for the bore; k = 0.5 and Poisson ratio 
o = 0.16974. 


IV. SURFACE-WAVE APPROXIMATIONS 


In this section, we consider two related equations that describe the 
high-frequency behavior of the surface-wave modes. We call these 
equations the lowest-order approximate equation and the refined ap- 
proximate equation. They are subject to the same restrictions about the 
disturbances being confined near the surface as are the asymptotic 
equations of Section III which describe the zeroth-order mode. The 
surface-wave approximations B.0 and B.1 (see Fig. 1) permit a more 
complete analysis of the higher-order modes. The refined approximation 
B.1 also describes the behavior of the modes in the transition region, at 
high frequencies, between the case of cross-sectional boundary curves 
of nonconstant (and not “almost” constant) curvature, for which the 
modes are localized, and the case of constant curvature, for which they 
are not localized. 

In the refined surface-wave approximation? the displacement, when 
evaluated at the surface of the cylinder, can be expressed in the form 
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Fig. 11—Asymptotic approximations 6p1, 5p2, and 6p3 to the normalized differential 
phase velocity as a function of the frequency parameter x for the wedge; k = 1.0 and Poisson 
ratio o = 0.16974. 
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(b2 + aa 1 dH 
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2b? G dp oy 
where H satisfies the refined approximate equation 
d*H 
as + {x[PK(n) — 2bv] — v2 + vSK(n) — 7[K(y)]?}H = 0. (16) 
U] 


Here the frequency parameter x is as defined in (2), K(n) is the curvature 
function, 7 = s/£ as before, and P, S, and 7 are constants. The parameter 
v is an eigenvalue, which is to be determined from a periodicity condition 
in the case of a closed boundary curve, and from an appropriate condition 
at infinity in the case of an open boundary curve. The propagation 
constant £ is given by 


Be = bx +», (17) 
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Fig. 12—Asymptotic approximations 6,2 and 5,3 to the normalized differential group 
velocity as a function of the frequency parameter x for the wedge; k = 1.0 and Poisson ratio 
o = 0.16974. 


and the normalized phase and group velocities are 
Wp = (b+ v/x)7}, Wg = (b+ dv/dx)7. (18) 


In this refined approximation, correction terms? of order 1/x could be 
included in the n and k components of the surface displacement given 
by (15). However, for the numerical cases considered in this paper, it 
turns out that these corrections, which differ for the two components, 
are of at most a few percent, so we do not write out these terms here. 

The lowest-order approximation B.0 (see Fig. 1) is obtained by 
omitting those terms multiplying H in (16) which are independent of 
x. Having done this, we replace H by Ho in (15) and v by vo in (17) and 
(18), where Ho satisfies the lowest-order approximate equation 





2H 
dn? + x[PK(n) — 2bv0]Ho = 0. (19) 


It was shown? that the asymptotic approximation (1) for the surface 
displacement of the zeroth-order mode may be derived from (15) and 
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Table l(a) — Asymptotic approximations to the normalized 
differential phase velocity of the zeroth-order mode for 
various cross-sectional shapes; frequency parameter 
x = 80 and Poisson ratio o = 0.16974 








Shape k —6p1 Xx 10° —dp2 XxX 108 —dp3 X 108 
Bore 0.3 —0.718 —1.072 —1.068 
0.4 —0.359 —0.767 —0.756 
0.5 0 —0.457 —0.442 
Rod 0.1 2.15 1.95 1.86 
0.2 2.51 2.22 2.10 
0.3 2.87 2.52 2.35 
0.4 3.23 2.82 2.61 
0.5 3.59 3.13 2.87 
Wedge 0.5 0.898 0.669 0.681 
1.0 1.80 1.47 1.43 
1.5 2.69 2.30 2.17 
Ridged 2 2.46 1.63 1.57 
Plane 3 3.69 2.67 2.45 
4 4.92 3.74 3.29 
5 6.15 4.83 4.09 





Table I(b) — Asymptotic approximations to the normalized 
differential group velocity of the zeroth-order mode for 
various cross-sectional shapes; frequency parameter 

x = 80 and Poisson ratio o = 0.16974 


Shape k —dy2 X 104 —6g3 X 104 
Bore 0.3 1.77 1.73 
0.4 2.04 1.92 
0.5 2.28 2.14 
Rod 0.1 1.02 1.83 
0.2 1.44 2.60 
0.3 1.77 3.32 
0.4 2.04 4.04 
0.5 2.28 4.78 
Wedge 0.5 1.14 1.02 
1.0 1.61 1.99 
1.5 1.98 3.18 
Ridged 2 4.17 4.65 
Plane 3 5.11 7.14 
4 5.90 10.10 
5 6.60 13.60 











the refined approximate equation (16). Also, the lowest-order asymptotic 
approximation, in which the terms involving C(n)[2(Px)1/2|-! do not 
appear in (1), may be derived from (19). However, in the transition region 
between the cases of nonconstant (and not “almost” constant) curvature 
and constant curvature, where the asymptotic results are not valid, eqs. 
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(16) and (19) have to be solved numerically, in general; it was shown,? 
though, that (19) may be solved analytically in the case of the curvature 
function K3(n) given by (6). We will see later how the results obtained 
from the numerical solution of (16) and (19) compare with the asymptotic 
results A.0 and A.1 (Fig. 1) in their common region of validity. 

As we will see, the lowest-order approximate eq. (19) is generally not 
sufficiently accurate in the transition region, when the curvature is 
“almost” constant. An exception is the case of the curvature function 
K3(n) given by (6), the reason being that the curvature is small in this 
case, tending to zero as k — 0, so that the terms involving S and 7 in (16) 
are small. 


4.1 Circular bore 


We first consider the case of a circular bore of radius @, with K(n) = 
—1, corresponding to k = 0 in (4). We compare our results from the 
lowest-order and refined approximate equations with the exact theo- 
retical results* for the fundamental mode in a circular bore. 

For the lowest-order surface mode, with K (ny) = —1 in (16) and (19), 
both Ho and H are constant, and the corresponding eigenvalues are 

2 1/2 
= — > —_ (ox +3) + | (ox +5) -@Px+7)| (20) 
Note that dvo/dx =0, so that in the lowest-order approximation, from 
(18), the group velocity is equal to the Rayleigh wave velocity. 

The exact theoretical dispersion relation for the fundamental mode 
in a circular bore was solved numerically for the normalized phase ve- 
locity w, by Rosenberg, Schmidt, and Coldren® for Poisson ratio 
o = 0.16974, corresponding to fused silica. They also calculated the 
corresponding value of the normalized group velocity wz. In Table II(a), 
we compare their values of w, as a function of the frequency parameter 
x with those calculated from (18) and (20). We add the subscripts L and 
R to w, to denote the lowest-order and refined approximate values of 
the phase velocity, respectively. Similarly, in Table II(b), we compare 
the exact theoretical and refined approximate values of w,. The nor- 
malized value of the Rayleigh wave velocity is 1/b = 0.905727. 

Notice that the lowest-order surface-wave approximation wpy, is 
reasonably close to w, and that the refined surface-wave approximations 
Wpr and Wgr are remarkably close to wp and wg, respectively. The 
agreement improves as the frequency parameter x increases. 

Rosenberg, Schmidt, and Coldren® plotted normalized differential 
phase and group velocities versus x. In Fig. 13, for purposes of compar- 
ison, the normalized quantities (bw —1)/(b — 1) are plotted against x 
for W = WpL, WpR, Wp, Wgr, and wy. The dots are for values corresponding 
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Table Il(a) — Comparison of lowest-order and refined approximate 
and exact theoretical values of the normalized phase velocity of 
the fundamental mode in a circular bore; frequency parameter 
x has various values and Poisson ratio o = 0.16974 








x WpL WpR Wp 

6 0.930315 0.938221 0.937969 

8 0.924044 0.928858 0.928829 
10 0.920321 0.923570 0.923583 
12 0.917856 0.920199 0.920219 
14 0.916104 0.917874 0.917893 
16 0.914794 0.916179 0.916195 
18 0.913777 0.914892 0.914905 
20 0.912966 0.913882 0.913893 
22 0.912303 0.913069 0.913078 





Table !I(b) — Comparison of refined approximate and exact 
theoretical values of the normalized group velocity of the 
fundamental mode in a circular bore; frequency 
parameter x has various values and 
Poisson ratio o = 0.16974 








x WsR Wy 

6 0.901046 0.902112 

8 0.902532 0.902825 
10 0.903408 0.903522 
12 0.903968 0.903962 
14 0.904347 0.904390 
16 0.904615 0.904644 
18 0.904812 0.904873 
20 0.904961 0.904967 
22 0.905077 0.905134 





to wp and wy. Note that the normalized values of wy do not lie precisely 
on asmooth curve. The values of w, were obtained through numerical 
differentiation once the values of w, had been calculated from the exact 
theoretical dispersion relation.® We suspect that the discrepancy is due 
to numerical difficulties in their computations. 


4.2 Ellipse-like bore and rod 


We now consider the ellipse-like bore and rod corresponding to the 
curvature functions K1(y) and K(n) given by (4) and (5). The eigenvalue 
problems for the refined and lowest-order approximate equations (16) 
and (19) were solved numerically. It suffices to consider the interval 
0 < n < 1/2, because the modes are either symmetric or antisymmetric 
about 7 =0, and about 7 = z/2, so that H’(0) = 0 or H(0) = 0, and 
H(x/2) = 0 or H(x/2) = 0. 
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Fig. 13—Normalized differential phase and group velocities (bw — 1)/(b — 1) asa 
function of the frequency parameter x for the circular bore, with Poisson ratio o = 0.16974. 
The curves correspond to the surface-wave approximations Wp, Wpr, and Wgr and the 
dots correspond to the exact theoretical results wp, and wy. 

A “shooting” method was used, which involves making an initial guess 
for the eigenvalue v, and numerically integrating the differential equation 
for H(n) from n = 7/2 to yn = 0. The value of v was adjusted iteratively, 
in the manner described in Appendix B, until the boundary condition 
at n = 0 was satisfied with sufficient accuracy. The initial iterations were 
done in single precision, and the final ones in double precision, and only 
a few iterations were required to obtain the desired accuracy. The nu- 
merical integrations were done from yn = z/2 to yn = 0, since in the as- 
ymptotic region the mode decays exponentially away from yn = 0, and 
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integration from 7 = 0 toward y = z/2 would lead to numerical insta- 
bilities in this case. 

To calculate both the phase and group velocities from (18), it is nec- 
essary to know the values of both v and dv/dy. Once the eigenvalue »p, 
and the corresponding eigenfunction H(n) had been obtained, the value 
of dv/dx was obtained by quadratures, using the expression (46) derived 
in Appendix B. The analogous expression for dvo/dx is given by (47). 


4.2.1 Rod 


We first compare the results of the numerical solution of (16) and (19) 
for the rod [k = 0.3 in (5)] with the asymptotic results A.0 and A.1 for 
x = 40. Here, and subsequently, the value of Poisson’s ratio is taken to 
be o = 0.16974. In Fig. 14, we plot the refined and lowest-order sur- 
face-wave approximations H(n) and Ho(n) for the zeroth-order sym- 
metric mode, normalized to unity at 7 = 0. The dots and circles corre- 
spond to the first and lowest-order asymptotic approximations S (7) 
and So(n), respectively, as defined in (10) and (9). It is seen that the as- 
ymptotic approximations agree quite well with the numerical solution, 
except, as expected, near n = 2/2. As the value of x increases, the 
agreement becomes better near 7 = 7/2, since the disturbance becomes 
exponentially small there. 

When the value of S,(7) is sufficiently small near n = 7/2, the 
zeroth-order mode which is antisymmetric about 7 = 7/2, as well as that 
mode which is symmetric about » = 7/2, is approximated by S,(m), as 
was argued in an earlier paper.” In Fig. 15, we plot H(m) for the lowest- 
order symmetric and antisymmetric modes. The dots, as before, corre- 





0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 
2n/7 


Fig. 14—Comparison of refined and lowest-order ee -wave approximations H(n) 
and H(n) for the zeroth-order symmetric mode, and corresponding asymptotic approxi- 
mations S;(n) (dots) and So(n) (circles) for the rod: k = 0.3, frequency parameter yx = 40, 
and Poisson ratio « = 0.16974 
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Fig. 15—Comparison of the refined surface-wave approximation H(y) for the lowest- 
order symmetric (S) and antisymmetric (A) modes, and the first-order asymptotic ap- 
proximation S;(7) (dots), for the rod; k = 0.3, frequency parameter x = 40, and Poisson 
ratio o = 0.16974. 


spond to the values of S;(n). To complete the comparison, in Table III 
we compare the asymptotic values of the differential phase and group 
velocities with the values obtained from (18). 

Trends to notice are that the asymptotic approximations 5,2 and 6,9 
to the differential phase and group velocities agree roughly with the 
lowest-order surface-wave approximations 6,; and 6,,. The corre- 
spondence between the next-order asymptotic approximations 6,3 and 
6,3 and the refined surface-wave approximations 6pr and dgr is somewhat 
better. The agreement is better for the differential phase velocities than 
it is for the differential group velocities. In all cases, the agreement im- 
proves as the frequency parameter y increases. 

As before, the t component of the displacement turns out to be a few 
percent of the n and k components of displacement. 

Higher-order modes may be investigated also by solving the eigenvalue 
problem (16) numerically. 


4.2.2 Bore 


The agreement between the asymptotic and the numerical results is 
even better for the bore. We have already discussed the circular bore and 
we now consider the transition from this to a noncircular bore for which 
the asymptotic results are good, by letting k vary from 0 to % in (4). The 
results of the numerical solution of the refined and lowest-order ap- 
proximate equations (16) and (19) for xy = 40 and o = 0.16974 are de- 
picted in Figs. 16 and 17 for the lowest-order symmetric and antisym- 
metric modes, respectively. The full curves give the values of H(n) (re- 
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Table Ill — Comparison of asymptotic and approximate values of 
the normalized differential phase and group velocities of the 
lowest-order symmetric (S) and antisymmetric (A) modes 
for the rod; k = 0.3, frequency parameter yx has 
various values, and Poisson ratio o = 0.16974 





























xX 

—6, X 103 40 60 80 
ose 5.745 3.830 2.873 
Spo 4.745 3.286 2.519 
dp (S) 4,783 3.3010 2.5271 

(A) 4.780 3.3007 2.5271 
Sng 4.089 2.9942 2.3550 
Spr (S) 4.075 2.9995 2.3555 

(A) 4.066 2.9919 2.3554 

—6, X 103 40 60 80 
by9 0.500 0.2729 0.1768 
du. (S) 0.428 0.2433 0.1611 

(A) 0.442 0.2448 0.1613 
bys 1.120 0.5475 0.3317 
bur (S) 1.188 0.5693 0.3401 

(A) 1.231 0.5730 0.3406 


fined approximation) and the broken curves the values of Ho(n) (low- 
est-order approximation); both curves are normalized to unity at n = 0 
for the specified values of k. 

For the symmetric mode, Fig. 16, Ho(n) = 1 and H(n) = 1 for k = 0, 
which agrees with the exact result for the circular bore. As k increases, 
the values of Ho(z/2) and H(z/2) decrease, becoming exponentially small 
for k = 0.5. It is seen that there is a significant difference between Ho(n) 
and H(n) for intermediate values of k. The lowest- and first-order as- 
ymptotic results agree very well with the numerical results obtained from 
(19) and (16) for k = 0.5. In Table IV we compare the lowest-order and 
refined approximations to the normalized differential phase and group 
velocities. We see that 6p, and dpr differ by only a few percent; the 
agreement improves as k increases. The lowest-order approximation 67, 
for the normalized differential group velocity, however, is not very good, 
especially for the smaller values of k. 

For the antisymmetric mode, Fig. 17, Ho(n) = cos 7 = H(n) for k = 0, 
which agrees with the exact result for the circular bore. The differences 
between Ho(n) and H(n), for intermediate values of k, are not as large 
as they are for the symmetric mode. It is noted that for k = 0.5, the curves 
of Ho(n) and H(n) are barely distinguishable from the corresponding ones 


98 THE BELL SYSTEM TECHNICAL JOURNAL, JANUARY 1977 


H(9), Hg (7) 





0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 
29/7 


Fig. 16—Refined and lowest-order surface-wave approximations H(y) (full curves) and 
Ho(n) (broken curves) for the lowest-order symmetric mode for the bore; k has various 
values, frequency parameter x = 40, and Poisson ratio o = 0.16974. 


for the symmetric mode, as expected from the asymptotic results. The 
most significant difference is that Ho(z/2) = 0 = H(z/2) for the an- 
tisymmetric mode, whereas Ho(x/2) = 0 = H’(z/2) for the symmetric 
mode. In Table V, we compare the lowest-order and refined approxi- 
mations to the normalized differential phase and group velocities. As 
was the case for the symmetric mode, 5,1 and bpp differ by a few percent 
and the agreement generally improves as k increases. The lowest-order 
approximation 6,7 for the normalized differential group velocity, al- 
though not very good quantitatively, is qualitatively better than for the 
symmetric mode. 
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Fig. 17—Refined and lowest-order surface-wave approximations H (yn) (full curves) and 


Ho(n) (broken curves) for the lowest-order antisymmetric mode for the bore; & has various 
values, frequency parameter x = 40, and Poisson ratio o = 0.16974. 
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Table IV — Comparison of lowest-order and refined 
approximations to the normalized differential phase and group 
velocities of the lowest-order symmetric mode for the bore; k 

has various values, frequency parameter x = 40, 


and Poisson ratio o = 0.16974 


k OpL OpR OgL OgR 
0 0.003605 0.003850 0 —0.224 X 1073 
0.01 0.003602 0.003847 —0.275 X 10-5 —0.226 X 1073 
0.02 0.003594 0.003835 —0.108 X 10-4 —0.234 X 1073 
0.03 0.003580 0.003817 —0.234 X 1074 —0.247 X 1073 
0.04 0.003562 0.003793 —0.396 X 10-4 —0.262 X 1078 
0.05 0.003539 0.003764 —0.583 X 10~4 —0.279 x 10-3 
0.1 0.003374 0.003560 —0.162 X 1073 —0.358 X 1078 
0.5 0.001232 0.001262 —0.581 X 1073 —0.636 X 1078 


Table V — Comparison of lowest-order and refined 

approximations to the normalized differential phase and group 

velocities of the lowest-order antisymmetric mode for the 
bore; k has various values, frequency 


parameter x = 40, and Poisson 


ratio 0 = 0.16974 


R OpL bp beh bgR 
0 0.003839 0.004066 —0.232 x 1073 —0.420 X 10-3 
0.01 0.003802 0.004024 —0.233 X 1073 —0.417 X 1078 
0.02 0.003764 0.003981 —0.235 X 1073 —0.415 X 1073 
0.03 0.003725 0.003936 —0.238 xX 1073 —0.414 x 1078 
0.04 0.003684 0.003890 —0.242 X 10-3 —0.415 X 1078 
0.05 0.003642 0.003842 —0.246 x 1073 —0.416 X 1078 
0.1 0.003419 0.003590 —0.279 X 1073 —0.430 X 1078 
0.5 0.001232 0.001262 —0.583 X 1073 —0.637 X 1073 
4.3 Wedge 


We now turn our attention to the wedge with a rounded tip, corre- 
sponding to the curvature function K3(n) given by (6). In this case, the 
lowest-order approximate equation (19) may be solved analytically.?:” 
The eigenfunctions are 


Ho(n) « (sech n)*"F [2am +m+1,-—m;am+1;%(1—tanhn)], (21) 
where 
Am = (Pxk + %4)1/2 -—-m—-—%>0,m=0,:--,M, (22) 
and the corresponding eigenfunctions are given by 
vo = a;,/(2bx). (23) 
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The hypergeometric function? in (21) terminates, and is a polynomial 
of degree m in its argument. There is a finite number of modes, because 
of the requirement that a,, > 0, which ensures that Ho(n) — 0 as 
|n| — «. The even-order modes are symmetric about 7 = 0, and the 
odd-order modes are antisymmetric, so it suffices to consider n = 0. 

According to (22), the zeroth-order mode, corresponding to m = 0, 
always exists for k > 0. For m = 0, the hypergeometric function in (21) 
is identically equal to 1. If Pyk < 2, in this approximation, then only the 
zeroth-order mode exists. In the limiting case, k — 0, corresponding to 
a planar boundary; ay — 0 and Hp — 1, for fixed yn. That is, the zeroth- 
order mode tends to a Rayleigh wave on a plane infinite half space as 
k—0. 

The eigenvalue problem for the refined approximate equation (16), 
with K(y) given by (6), was solved numerically by a shooting method 
after a transformation of the independent variable had been made to 
reduce the interval of integration to a finite one. The details are given 
in Appendix C. For the symmetric modes, H’(0) = 0, and for the an- 
tisymmetric modes, H(0) = 0. Once the eigenvalue v, and the corre- 
sponding eigenfunction H(n) had been obtained, the value of dv/dy was 
obtained by quadratures, using the expression (53) derived in Appendix 
C. The values of v and dv/dx were used in (18) to obtain the normalized 
phase and group velocities. 

In Fig. 18, we plot the refined and lowest-order surface-wave ap- 
proximations H(n) and Ho(n), normalized to unity at 7 = 0, for the 
zeroth-order mode for the wedge for k = 1, x = 40, and o = 0.16974. We 
also make a comparison with the asymptotic results, A.0 and A.1 (Fig. 
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Fig. 18—Comparison of refined and lowest-order surface-wave approximations H(n) 
and Ho(n) for the zeroth-order mode, and corresponding asymptotic approximations S1(y) 
(dots) and So(n) (circles) for the wedge; k = 1.0, frequency parameter x = 40, and Poisson 
ratio o = 0.16974. 
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1), the dots and circles corresponding to the first- and lowest-order as- 
ymptotic approximations S,(7) and So(), respectively, as defined by 
(10) and (9). In Table VI, we compare the asymptotic values of the dif- 
ferential phase and group velocities with the values obtained from (18). 
The trends are similar to those we have noticed for the other cases. The 
asymptotic approximations 6p and 62 agree roughly with the lowest- 
order surface-wave approximations 6), and 46g,. The higher-order 
asymptotic approximations 6,3 and 6,3 agree better with the refined 
surface-wave approximations 6pr and égr. The convergence is better for 
the differential phase velocities than it is for the differential group ve- 
locities; also, the agreement between the asymptotic and the surface- 
wave approximations is better. 

In Figs. 19 and 20 we plot H(n) and Ho(n) for the remaining three 
modes. The odd-order, antisymmetric modes are normalized so that 
H’(0) = 1 and Ho(0) = 1. In Table VII, we compare the corresponding 
values of the differential phase and group velocities. 

Finally, we consider the transition region, between the case of non- 
constant (and not “almost” constant) curvature and constant curvature 
for the wedge, with x = 40 and o = 0.16974. In Fig. 21, we plot H(n) for 
the zeroth-order mode for several values of k between 0.01 and 1. We 
have not plotted Ho(n) in this figure, since we compared Ho(n) with H(n) 
in Fig. 18 for k = 1 and found differences to be quite small for the smaller 
values of k. This is because the curvature is small when k is small, tending 
to zero as k — 0, so that the terms involving S and 7 in (16) are small. In 
Table VIII, we compare the lowest and refined approximations to the 
differential phase and group velocities as obtained from (18). In the 
transition region the agreement is good, both for the differential phase 
velocity and for the differential group velocity. As k increases and the 
lowest-order approximate equation becomes less accurate, discrepancies 
appear. 


Table VI — Comparison of asymptotic and approximate values of 
the normalized differential phase and group velocities of the 
zeroth-order mode for the wedge; k = 1.0, frequency 
parameter x = 40, and Poisson 
ratio 0 = 0.16974 








—5, X 103 x = 40 —5, X 103 x = 40 
5p1 3.591 = = 
Spo 2.678 by 0.457 
SL 2.778 Set 0.351 
dps 2.515 by 0.605 
5pR 2.521 bk 0.614 
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Fig. 19—Refined and lowest-order surface-wave approximations H(n) and Ho(n) for 
the second- aaa mode for the wedge; k = 1.0, frequency parameter xy = 40, and Poisson 
ratio o = 0.16974 


V. SUMMARY 


We made numerical computations to learn about the propagation of 
elastic surface waves along cylindrical objects roughly corresponding 
to an elliptical bore, an elliptical rod, a wedge with a rounded tip, and 
a flat plane with a rounded ridge. The cross-sectional curvature functions 
describing these objects are given by eas. (4) to (7). In earlier papers,” 
we had derived two approximate analytical descriptions of the surface- 
wave behavior: a high-frequency asymptotic approximation A, and one 
that we termed a surface-wave approximation B, as depicted in Fig. 1. 
Each of these approximations, in turn, was available in two forms: a 
lowest-order one and a higher-order or refined one. Here, we evaluated 
these approximations numerically. 

We first performed a high-frequency asymptotic analysis of the dis- 
turbance in the vicinity of the cylinder surface and obtained the low- 
est-order, A.0, and next-higher-order, A.1, asymptotic approximations.” 
We used these approximations in the form shown in eq. (1), which de- 
scribes the zeroth-order mode at the surface of the cylinder. For the bore 
and the rod, this equation corresponds to both the zeroth-order sym- 
metric and antisymmetric modes. We also used the high-frequency 
asymptotic approximations to the phase and group velocities given by 
(12) and (13). The analysis involved two restrictions: the frequency had 
to be high enough that the disturbance was confined close to the surface 
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Fig. 20—Refined and lowest-order surface-wave approximations H(y) and Ho(n) for 
the first- and third-order modes for the wedge; k = 1.0, frequency parameter x = 40, and 
Poisson ratio o = 0.16974. 
of the cylinder, and the deviation of the cross-sectional curvature from 
a constant value had to be sufficiently large that the disturbance was 
confined near points of maximum algebraic curvature. 

We also were able to describe the mode behavior at the cylinder sur- 
face? by the lowest-order, B.0, and refined, B.1, surface wave approxi- 
mations (16) and (19). The phase and group velocities were given by (18). 
These surface-wave approximations B were subject to the same fre- 
quency restriction as were the asymptotic approximations A. In fact, the 
lowest-order and next-highest-order asymptotic approximations A.0 
and A.1 for the zeroth-order mode could be obtained from the lowest- 
order and refined surface-wave approximations B.0 and B.1, respectively, 


Table Vil — Comparison of the lowest-order and refined 
approximations to the normalized differential phase and 
group velocities of the four modes for the wedge; 

k = 1.0, frequency parameter x = 40, 
and Poisson ratio o = 0.16974 


m —dpL xX 10° —OpR xX 10° —OoOgL X 103 —OgR x 103 
0 2.778 2.521 0.351 0.614 
1 1.408 1.231 0.821 0.944 
2 0.497 0.394 0.828 0.836 
3 0.050 0.021 0.370 0.264 
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H (7) 





Fig. 21—Refined surface-wave approximation H(n) for the zeroth-order mode for the 
wedge; various values of k, frequency parameter x = 40, and Poisson ratio o = 0.16974. 


provided that the above restriction on the cross-sectional curvature was 
satisfied. 

The advantage of the surface-wave analysis B was that the curvature 
restriction could be dropped. It was possible, with the refined surface- 
wave approximation B.1 to describe the mode behavior, at high 
frequencies, in the transition region between the case of cross-sectional 
boundary curves of nonconstant (and not “almost” constant) curvature, 
for which the modes are localized, and the case of constant curvature, 
for which they are not localized. Also, the surface-wave approximations 
B permitted a more complete analysis of the higher-order modes. A 
disadvantage of the surface-wave approximations B was that they con- 


Table VIIl — Comparison of lowest-order and refined 
approximations to the normalized differential phase 
and group velocities of the zeroth-order mode 
for the wedge; k has various values, 
frequency parameter x = 40, and 
Poisson ratio o = 0.16974 





k —dpi X 103 —dpr X 108 —dy, X 104 —dyp X 10! 
0.01 0.04302 0.04294 0.03381 0.03384 
0.02 0.1426 0.1422 0.09537 0.09561 
0.03 0.2763 0.2750 0.1635 0.1643 
0.04 0.4333 0.4308 0.2324 0.2341 
0.05 0.6076 0.6034 0.3003 0.3033 
0.1 0.1639 0.1620 0.6116 0.6278 
0.25 0.5424 0.5284 1.337 1.466 
0.5 1.254 1.194 2.291 2.801 
1.0 2.778 2.521 3.513 6.137 
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sisted of eigenvalue equations which, in general, had to be solved nu- 
merically. For the particular boundary curves chosen, the asymptotic 
approximations A involved little more than numerical evaluation of some 
analytical formulas, once the quadratures had been done analytically. 

In general, we found numerically that the lowest- and next-higher- 
order asymptotic approximations A.0 and A.1 did agree with the low- 
est-order and refined surface-wave approximations B.0 and B.1, re- 
spectively, in their common region of validity. This was true both of the 
results for the disturbances, and for the phase and group velocities. The 
agreement improved as the frequency parameter x was increased, and, 
for the phase and group velocities, was better between the higher-order 
and refined approximations than it was between the other two. 
The asymptotic and surface-wave approximations for the disturbance 
did not agree particularly well for the case of cylinders with closed 
boundary curves for values of 7 for which the disturbance was expo- 
nentially small. This was to be expected, since one expression had to 
suffice for both the lowest-order symmetric and antisymmetric modes 
in the asymptotic approximation A, while separate expressions were 
available in the surface-wave approximation B. We also used the refined 
surface-wave approximation B.1 numerically to describe disturbances 
in the transition region discussed earlier, and used the lowest-order and 
refined surface-wave approximations B.0 and B.1 to investigate the 
higher-order modes. 

We turn now to a qualitative description of the numerical results. We 
first discuss the phenomenon of mode confinement and its dependence 
upon such things as the shape of the cylinder and the value of the fre- 
quency parameter x. We then discuss our results for the phase and group 
velocities. 

We found that the t component of the displacement for the lowest- 
order mode was only a few percent of the size of the n and k components. 
These latter two, when normalized to unity at n = 0, were either the same 
as a function of n (asymptotic theory A) or differed by a few percent at 
most (surface-wave theory B). It thus sufficed to consider a normalized 
scalar displacement function, rather than a vector function. The com- 
plete solution is essentially that for Rayleigh waves traveling on the 
surface of a plane infinite half space except that it is multiplied by a 
function of 7, which describes the confinement of the wave due to the 
cylinder curvature, or, to be more precise, the confinement due to the 
deviation of the cylinder curvature function from a constant value. It 
is this confinement function (the normalized scalar displacement 
function) that we computed. 

The wedge with a rounded tip and the plane with a ridge on it are 
cylinders whose cross-sectional boundaries are each described by a 
curvature function with a single algebraic maximum at 7 = 0. For these 
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cylinders, we found that the confinement function decays rapidly with 
|n| away from its value of unity at 7 = 0. This means that the surface- 
wave disturbance is confined to the vicinity of the tip of the wedge or the 
ridge on the plane. In both cases, the amount of confinement decreases 
as the parameter k decreases and the cylinder becomes more nearly 
planar. For the ridged plane, we observed this when the curvature 
function was not “almost” constant (i.e., for k not too small). For the 
wedge, we also made computations for small values of k, corresponding 
to the flattening out of the wedge into a plane. Here we were aided by 
an analytical solution of the lowest-order surface wave approximation 
B.0. It showed that, as k — 0, only the zeroth-order mode exists, that it 
is symmetric, and that it tends to a Rayleigh wave on a plane infinite half 
space. This was confirmed by the numerical computations, both in the 
lowest-order and in the refined surface-wave approximations B.0 and 
B.1. 

We also plotted the confinement functions for higher-order modes 
on the wedge. For given values of k and the frequency parameter y, there 
are finitely many surface modes. The even-order modes are symmetric 
about n = 0, and the odd-order modes are antisymmetric. 

The ellipse-like rod and bore have boundary curves that are symmetric 
and which attain their algebraic maxima at two points, n = 0 and n =z. 
We investigated only the zeroth-order modes, although higher-order 
modes may also be studied numerically. The asymptotic approximation 
A for the zeroth-order mode on a rod or bore actually corresponds to two 
modes, a symmetric one and an antisymmetric one. These can be treated 
separately with the surface-wave approximation B. For a cylinder whose 
curvature is not “almost” constant, we observed confinement of the 
displacement to two regions. Each cylinder cross-section has two points 
of maximum algebraic curvature. They define two generators of the 
cylinder. The displacement is confined in the vicinity of these genera- 
tors. 

We considered the transition from an ellipse-like bore with the definite 
confinement properties discussed above to a circular bore (k = 0), which 
exhibits no confinement at all. For small values of k, it was necessary to 
use the refined surface-wave approximation B.1 rather than the low- 
est-order approximation B.0 in order to describe the modes adequately. 
We treated the lowest-order symmetric and antisymmetric modes. For 
k = 0, the results agreed with the known analytical results for a circular 
bore: the confinement function is constant for the symmetric mode and 
goes like cos 7 for the antisymmetric mode. As k increased, confinement 
began to appear. As k approached 0.5, there was definite confinement 
and the surface-wave approximate results B agreed with the asymptotic 
results A, which had been obtained earlier and which were valid for 
curvature functions that were not “almost” constant. 
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These were the basic confinement properties that we observed when 
we studied cylinders described by several different cross-sectional cur- 
vature functions with a variable parameter k. We also varied the fre- 
quency parameter x and found that, for any given cylinder, the con- 
finement becomes more pronounced as x increases. 

For most of the numerical computations we chose the Poisson ratio 
to be « = 0.16974, corresponding to fused silica. A few computations were 
made with other values of o; we found that the confinement increases 
as o decreases. 

We also calculated the phase and group velocities. In the asymptotic 
approximation A, these are given by explicit asymptotic formulas. In 
the surface-wave approximation B, the velocities are given in terms of 
an eigenvalue and its derivative with respect to x. The eigenvalue was 
determined from a periodicity condition in the case of a closed boundary 
curve and from an appropriate condition at infinity in the case of an open 
boundary curve. The derivative of the eigenvalue with respect to x was 
expressed in terms of quadratures, which were evaluated numerically. 
This avoided the difficulty of numerical differentiation with respect to 
the frequency. Both the phase and group velocities tend to the Rayleigh 
wave velocity as x > ~. We computed the differential phase and group 
velocities normalized with respect to the transverse-wave velocity cr. 

The trends that we generally observed were that the asymptotic ap- 
proximations 6p2 and 6,2 to the differential phase and group velocities 
agreed roughly with the lowest-order surface-wave approximations dpz, 
and 6,,. Better agreement was obtained between the next-order 
asymptotic approximations 5,3 and 6,3 and the refined surface-wave 
approximations 6pr and dsr. The convergence and the agreement im- 
prove as the frequency parameter x increases. 

In the transition region between cylinders of constant curvature and 
those of not “almost” constant curvature (where the parameter k is small 
and the asymptotic theory is not valid), the lowest-order surface-wave 
approximation B.0, as expected, was not always too good, particularly 
for the differential group velocity, so it is necessary to use the refined 
surface-wave approximation B.1. 

Finally, we compared the surface-wave approximation results B for 
the circular bore with exact theoretical results and obtained excellent 
agreement. 
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APPENDIX A 

We summarize here the asymptotic results? in terms of suitably nor- 
malized quantities. The longitudinal and transverse velocities are given 
by cy = [(A + 2u)/p]!/2 and ep = (u/p)1/2, where \ and u are Lamé’s con- 
stants. Also, the Poisson ratio is o = \[2(\ + u)|~!. The normalized 
Rayleigh wave velocity is wr = cr/cr, where wp is the root of the equa- 
tion 


(24) 


(1 — uy = 1 — wip [1-5 wit |, 


ia) 


which satisfies 0 < wr < 1. We define the quantities 


= (1 — 2c) 


b=1 =| b2 
/wp, a | d(1. = ¢) 


1/2 
| ape@zsaye. 0b) 


and 
2 b2 = : 
Pee (26) 
b*(ay, — ar)? + 2aza7(ay — ar) 
We further define the quantities R and 7 by means of the equations 


P2(a, —a 


2 
4(a, — ar) [b*(a, — ar) + 2a,a7]R = Da? r) [b2(a, + ar)? 
LYT 


— dajaz]+ 2Pl[az(b?2 — azar) — ar(az — ar)?| 
+ [b2(a% — 8a7) + 2azaz], (27) 


and 


(a, —ar)[b%(ay — ar) + 2aza7z](7 + R) 





b2 
=f | (ai + a} — azar) + azar(2azy — Sar) | (28) 
atatT 


If » is small, the curvature function K() has an expansion of the 
form 


K(n) = do + don? + dgn? + dant +---. (29) 


There is no term proportional to 7 in (29), because of our assumption that 
the curvature attains its algebraic maximum at n = 0. It is assumed that 
dy <0. We define the quantity 


11 sd3\2_ 3d4_— /doPy? 
=— (2) —-—*_ (=) + am. 30 
Origa) cada) ” 
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Then, from the asymptotic results,? the reciprocals of the normalized 
phase and group velocities defined in (11) have the expansions 


14 dP (doP2 Q 








. ve 31 
Wp 2bx 2bx3/2 2bx2 (31) 
and 
1 (-doP)'?_—Q 
ie pepe oa 32 
Wg Aby 3/2 2byx? 2) 


If we expand the reciprocals of these expansions we obtain those given 
in (12) and (18). 

We now define the functions C(y), F(n), G(n), and J(n) occurring in 
(1). In terms of the curvature function K(n), we define 


I(n) = [do — K(y)]}” sgn 7, (33) 
and let 


[K’(n) + 2(—ds)/2I(m)] 


Tay = 
2 4[do — K(n)] 


(34) 


The prime denotes differentiation with respect to the argument, and it 
is seen from (29) that L(0) is finite. Then, we define 


Oe | f "Leds, Gay= f" 16ae. (35) 


Next, we let 


M(n) = f, "TOK(Ods, (36) 


and 
N(n) = §, MIL) + (LIQ)? — L'(0) — [LO)PVIOdk 87) 
Finally, we define 
C(n) = N(n) — doRG(n) + 7M(n). (38) 


Then, from the asymptotic results,” the disturbance corresponding to 
the lowest-order mode can be expressed as 
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cet (bz—wt) 


wh (0) 
2 2 

= F(n) exp [-(Px)!2G(n)] (see 

T 


24 92 
x E + Tone n- | evouxe erent | 


u(Z,n) 


-aTX= — apeneix | 


2(Px)1/2 2b2 
«{(—) "Ton + ib E rs A k|). (39) 


APPENDIX B 


We describe here the “shooting” method used to solve numerically 
the eigenvalue problems for the refined and lowest-order approximate 
equations (16) and (19) for the bore and rod, corresponding to the cur- 
vature functions K (ny) and Ko(7) given by (4) and (5). It is desirable to 
shoot from n = z/2 to n = 0, since in the asymptotic region the mode 
decays exponentially away from n = 0, and integration from y =-0 toward 
n = x/2 would lead to numerical instabilities. Consequently, we let 


E = 1/2 —, Z\(E) = H(y), Zolé) = dZ,/dé. (40) 


Since the eigenvalue v has to be determined, we also consider the dif- 
ferential equations for 


OZ 
ov’ 
The initial conditions are taken as 


Z,(0) = 1, Z(0) = 0, or Z;(0) = 0, Z(0)=1, (42) 


AC ee ee (41) 
Ov 


according to whether the mode is symmetric, or antisymmetric, about 
n = 1/2. In either case, the remaining initial conditions are 


Z3(0) = 0, Z4(0) = 0. (43) 


An initial guess for the value of vy was made, and the system of equations 
for Z;(£), i = 1,2,3,4, was integrated from £ = 0 to = 2/2. When a mode 
symmetric about 7 = 0 was sought corresponding to Z2(a/2) = 0, the 
initial value of »y was changed to v — Zo(a/2)/Z4(a/2), since Zo(a/2, vp + 
6) © Zo(x/2,v) + 60Z2/dv(7/2,v). Analogously, if a mode antisymmetric 
about 7 = 0 was sought, corresponding to Z,(a/2) = 0, then the initial 
value of v was changed to v — Z1(a/2)/Z3(a/2). The process of integrating 
the system of equations was repeated, until the boundary condition at 
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n = 0 was satisfied with sufficient accuracy. The initial iterations were 
done in single precision, and the final ones in double precision, and only 
a few iterations were required to obtain the desired accuracy. 

The shooting program was checked in the case of the lowest-order 
approximate equation (19), since the eigenfunctions may be expressed 
in terms of Mathieu functions? when the curvature function is given by 
(4) or (5). The checks were carried out for values of the parameter g = 
Pxk equal to 1, 5, and 10, the eigenvalues and eigenfunctions being 
checked against tabulated values.!° 

We now turn our attention to the calculation of dv/dy, which is needed 
to calculate the group velocity from (18). If we let H, = 0H/dx, then from 
(16) we obtain 
ay + {y[PK b 2+ ySK K(n)|33H 

dn? {x [PK (n) — 2bv] — v? + vSK(n) — 7[K(n) |} x 


= {[2(bx + v) — SK(n)] a + [2bv — PK(n)]} H(n). (44) 


Hence, 

d dH dH d?H. d?H 
Beit 2 gk eam 8 a) ys peel) ay 2 genet! 
dn ( dn * dn dn? * dn? 


= [20x +v) —SK(n)] o + [2b — PK(»)]| [H(n)]?. (45) 


But from the boundary conditions at n = 0 and n = 2/2, it follows that 
[H dH,/dn — H,, dH/dy]s’” = 0. Hence, if we integrate (45) from 7 = 0 
to n = 1/2, we obtain 

dv 


n/2 
dx J, [2(bx + v) — SK(n)|[H(n)]?dn 


a/ 
-f * [PK(n) — 2bv|[H(m)2dn. (46) 


Analogously, from (19), it follows that 
dvo 

‘s dx Jo 

The program was written so that the system of equations for Z;(&), 1 = 


1,2,3,4, was augmented in the double-precision stage of the iterations 
to include the evaluation of the two quadratures in (46) or (47). 


x/2 x/2 
2b [Ho(n)|2dn = f, [PK (n)— 2bvollHo(n)|2dn. (47) 


APPENDIX C 


We describe here the “shooting” method, similar to that described 
in Appendix B, that was used to solve numerically the eigenvalue 
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problem for the refined approximate equation (16) for the wedge, cor- 
responding to the curvature function K3(7n) given by (6). From symmetry, 
it suffices to consider the interval n = 0. It is desirable to make a trans- 
formation of variables that reduces this interval to a finite one, partic- 
ularly since we want to integrate from 7 = ~ toward 7 = 0, in order to 
avoid numerical instabilities in the asymptotic region. Consequently, 
we introduce the new independent variable 


£= 'b(1 — tanh y), (48) 


which is suggested by the form of the solution (21) of the lowest-order 
approximate equation (19). This form also suggests the substitution 


H(n) = (sech n)*g(§), aw = (2bxv + v?)¥?2 > 0. (49) 
From (6), (16), (48), and (49), it follows that 


_ Ws op 8 
KIS) ari-20 


+ [(Px + vS)k — a(a + 1) — 4rk2&(1 — Olg =0. (50) 


The range of fis from 0 to }4, with ¢ = 0 corresponding to n = ©. Exam- 
ination of the behavior of the solutions of (50) for ¢ — 0, and the re- 
quirement that H(n) — 0 as 7 > ~, lead to the condition that g(0) be 
finite. The value of g’(0) may be determined by setting ¢ = 0 in (50). 
Thus, as initial conditions, we take 


g(0)=1, (0) =a — (Px + vS)R/(a + 1). (51) 


Since the eigenvalue v has to be determined, we also consider the dif- 
ferential equation for 0g/dv. 
Because the coefficient of d2g/d£? in (50) vanishes at ¢ = 0, we let 


Yi =g-1, Yo = dY,/d¢é — g’(0), Y3 = OY;/0d», 
YY, = OY>2/d», (52) 


so that Y;(0) = 0, and Y;(0)/¢ is finite at ¢ = 0, i = 1,2,3,4. The system of 
equations for Y;({) was integrated from ¢ = 0 to [= }4, the value of v being 
adjusted after each step of the iteration procedure until the condition 
g’(44) = 0, or g(4) = 0, was satisfied with sufficient accuracy. The former 
condition corresponds to a mode that is symmetric about n = 0, and the 
latter to one that is antisymmetric. 

It remains to discuss the calculation of dv/dx. If we integrate equation 
(45) from n = 0 to y = &, with K(n) given by (6), it follows that 
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d co 
rh f, [2(by + v) — Sk sech? n][H(n)|2dn 


= f ” (Pk sech? n — 2bv)[H(n)|2dn. (53) 


In terms of the new variables given by (48), (49), and (52), this requires 
the evaluation of the definite integrals 


/ 
f° asa - diet + YP, 


/2 
S [45 — Die V(O[2 + Vids, (64) 
and the calculation of 


/2 
f {4¢(1 a Oletdé= _VaT (a) 


The integral in (55) was expressed in terms of gamma functions,!! for 
which a double-precision routine was available, in order to avoid a sin- 
gular integrand at ¢ = 0 when 0 < a < 1. The integrals in (54) were 
evaluated in the double-precision stage of the iterations by augmenting 
the system of equations for Y;({), i = 1,2,3,4. 

To check the accuracy of the shooting method, the analogous system 
of equations corresponding to the lowest-order approximate equation 
(19) was solved numerically, and the results were checked against those 
calculated from the analytical solution (21) to (28). 
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